CC BY
14
УДК 550.34: 004.41
АЛГОРИТМ И ПРОГРАММЫ МОДЕЛИРОВАНИЯ ДВУХМЕРНЫХ СЕЙСМИЧЕСКИХ ПОЛЕЙ И ИХ ПРАКТИЧЕСКОЕ ПРИМЕНЕНИЕ
Дмитрий Алексеевич Караваев
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Лаврентьева, 6, кандидат физико-математических наук, научный сотрудник, тел. (383)330-70-69, е-mail: kda@opg.sscc.ru
Представлен алгоритм численного моделирования распространения упругих волн в двухмерных неоднородных изотропных упругих средах. Разработана и предложена параллельная реализация для возможности проведения расчетов на высокопроизводительных вычислительных системах с различной архитектурой. Разработаны программы для проведения расчетов на основе конечно-разностного метода четвертого порядка аппроксимации по пространству. Представленная технология математического моделирования с использованием MPI и CUDA применена для расчета реалистичной модели по сейсмическому профилю на 400 км. Представлены синтетические данные о структуре сейсмического поля для Байкальской рифтовой зоны.
Ключевые слова: параллельный алгоритм, программа, сейсмические волны, большие расстояния, численное моделирование, реалистичные модели.
ALGORITHM AND PROGRAMM CODES FOR 2D SEISMIC FIELD SIMULATION AND COMPUTATION, AND THEIR PRACTICAL APPLICATION
Dmitry A. Karavaev
Institute of the Computational Mathematics and Mathematical Geophysics SB RAS, 6, Prospect Аkademik Lavrentiev St., Novosibirsk, 630090, Russia, Ph. D., Researcher, phone: (383)330-70-69, е-mail: kda@opg.sscc.ru
An algorithm for simulation of elastic wave propagation in 2D inhomogeneous isotropic elastic media is presented. The parallel realization of algorithm to perform computations on high-performance computing systems with different architectures is developed and proposed. Program codes to perform calculations on the basis of finite-difference method of fourth order respect to space were developed. The presented technology of mathematical modeling with usage of MPI and CUDA is used to perform computations for long distance realistic model for seismic profile of 400 km of size. The obtained new synthetic data depicting the structure of seismic filed for the Baikal rift zone are presented.
Key words: parallel algorithm, program code, seismic waves, numerical simulation, realistic models.
Введение
Теоретические и экспериментальные исследования распространения сейсмических волн в неоднородных средах со сложной геометрией строения на большие расстояния являются весьма актуальными [1-3]. В работе основным объектом для изучения выступает модель Байкальской рифтовой зоны [4]. Такой объект интересен тем, что известны несколько геофизических моделей, какими их представляю геологи. Проведение математического моделирования распространения упругих волн в такой неоднородной среде позволит получить
новые знания о его строении и возможно поможет уточнить строение модели. Такое исследование возможно при проведении сравнения экспериментальных и синтетических результатов.
Разработка алгоритма, его параллельной реализации и программного обеспечения вызвана необходимостью проведения расчетов по профилю на расстояния сотни километров. В этой связи важными задачами выступают поиск и разработка математического метода с порядком точности, обеспечивающим приемлемую точность вычислений. В настоящее время существуют разнообразие математических методов, которые применимы для решения динамической задачи теории упругости. Это могут быть интегральные, конечно-разностные методы, методы конечных элементов [5-9].
Современные задачи в области численного моделирования, связанные с расчётами и сеточными методами как правило являются весьма трудоемкими. В таких задачах исследователь сталкивается с необходимостью работы с большими объемами данных. В таких условиях становиться необходимым разработка параллельных алгоритмов и применение для расчетов высокопроизводительных многоядерных вычислительных систем различной архитектуры [10-15]. Такие системы при правильном и эффективном использовании вычислительных ресурсов [16-18] на основе параллельных алгоритмов позволяют проводить исследования для реалистичных моделей с подробной детализацией. Необходимость проведения расчетов на кластерах объясняется тем, что речь идет о работе с десятками и сотнями Гб данных.
В работе представлена технология проведения математического моделирования для расчета полного волнового поля для неоднородных упругих сред. Выбран конечно-разностный метод четвертого порядка аппроксимации. Для исключения отражений волн от границ расчетной области применен алгоритм PML (Perfectly Matched Layers) [19-24]. Для выбранного метода предложено распараллеливание численного алгоритма. На основе параллельного алгоритма разработано программное обеспечение для кластеров с различной архитектурой. Речь идет об использовании в расчетах центральных процессоров, графических процессоров. Предложенная параллельная реализация легко адаптируется к использованию различных вычислительных устройств. На основе разработанных программ проведена работа по восстановлению структуры сеточной модели и определения значений упругих параметров для проведения численного моделирования сейсмического поля рифтовой зоны. В данной работе впервые получены и представлены результаты моделирования полного волнового поля для Байкальской рифтовой зоны, которые не имеют аналогов.
Постановка задачи
В работе рассматривается решение прямой задачи геофизики численного моделирования полного волнового поля для изотропного случая неоднородных упругих сред. Для решения задачи используется система уравнений теории упругости (1), записанная в терминах скоростей перемещений u, w и напряжений
тхх ,т22 ,гХ2. Решение ищется в двухмерной прямоугольной системе координат (х, 2). Область моделирования представляет собой прямоугольник, верхняя граница которого (2=0) является свободной поверхностью, рис. 1.
/ свободная поверхность
РМ1.
РМЬ
Рис. 1. Схематическое представление 2Б области моделирования
Упругая среда характеризуется тремя параметрами: плотностью р, скоростями продольных и поперечных волн Ур и Ув. Геометрия строения внутренней части прямоугольной области определяется распределением значений представленных упругих характеристик. Отметим, что в системе уравнений теории упругости эти характеристики входят как параметры Ламе. Эти параметры являются двухмерными функциями координат и определяются пользователем до начала расчетов волнового поля.
Риг = Тхх, Х + ТХ2,2 ,
Р^г =*Х2, Х ,2 ,
Тхх, = (Я + 2М)иХ + МЩ + /х , (1)
Т22 ,г =Мих + + 2МК + Л ,
^ ,г = Киг + х )
= 5м = 3^ = =
где и( =— ,Ту,. =— , Х1 = х, Х2 = 2
дг дх]
Задача решается при соответствующих нулевых граничных и начальных условиях.
Частью успешного решения является применение метода РМЬ в одной из его реализаций для исключения нежелательных отражений упругих волн от границ расчетной области, которые могли бы внести изменения в интересующую картину сейсмического поля внутри области моделирования [25, 26]. Для этого вдоль боковых границ и вдоль нижней границы выделены прямо-
угольные подобласти малого размера, в которых и реализуются расчетные формулы метода PML, рис. 1.
Для проведения расчетов и генерации сейсмических волн используется точечный источник, который может быть реализован как «всесторонняя сила», находясь внутри области моделирования исключая PML слои, или как «вертикальная сила», располагаясь на свободной поверхности.
Метод решения
При выборе метода решения было важно, чтобы он обладал хорошей гибкостью и потенциалом для разработки параллельного алгоритма. Также метод должен обеспечивать достаточную точность вычислений на больших расстояниях. В этой связи был выбран конечно-разностный метод четвертого порядка точности по пространству [27] с уравнениями для реализации граничных условий на свободной поверхности [28]. Конечно-разностная схема для расчета P-SV волн имеет общий вид (2):
С1,, = С,,, + £ (К + 2м\J [ А(и
3
1 +-, j 2
и
1 - ^
) + В(ип 1 - ип 1 )] +
1+-, J 2
1-^
А [[А(- 3 - -п 3 ) + В(-п 1 - -п! ) + /хп+1]), А = -—,В = -
1
(2)
1+-, J 2
1- ^
1+-, J 2
1-^
24
8
где 1, j - сеточные координаты, п, п +1- временные слои, At, Ах - шаг расчетной сетки по времени и пространству соответственно.
Схема реализована на смещенных сетках [28], что отражается на своеобразном размещении компонент вектора скоростей смещений и компонент тензора напряжений относительно выбранного узла сеточной модели, рис. 2.
Рис. 2. Схема на смещенных сетках для 2D случая. и,— - скорости перемещений, тхх,тгг,тхг - напряжения
Для данной работы расположение выбрано таким образом, что на свободной поверхности находятся компоненты txx,tzz,u. В работе рассматриваются сеточные модели только с одинаковыми шагами по пространственным переменным Ax = Az.
Программная реализация
Для данной работы разработаны несколько программных реализаций на основе параллельного алгоритма. Для проведения расчетов применяются кластеры на основе многоядерных CPU (Central Processing Unit) и на основе специализированных вычислительных устройств типа GPU (Graphics Processing Unit). Выбран простейший способ организации параллельных вычислений на основе разбиения по данным. Таким образом область моделирования представляется как набор горизонтальных слоев, выбран способ одномерной декомпозиции по оси Oz. Каждый слой обрабатывается независимо на выделенном вычислительном устройстве. Разработанная параллельная реализация алгоритма легко адаптируется для использования в расчетах CPU или GPU. Таким образом программа состоит из нескольких частей. Первая часть реализует организацию параллельных вычислений, выделение процессов, создание топологии, в данном случае «линейка», интерфейс взаимодействия между вычислительными устройствами, организацию пересылок данных, ожидание выполнения операций и их проверку. Все это реализовано с использованием технологии MPI (Message Passing Interface). Вторая часть программы реализует процедуры расчета, которые могут быть распараллелены с использованием OpenMP (Open Multi-Processing) или технологии CUDA (Compute Unified Device Architecture) [29]. Поскольку CUDA в основном оперирует с одномерными массивами данных при реализации распараллеливания, то в разработанных программах было реализована конвертация двухмерных массивов в одномерные. Также при разработке программ следует учитывать возможности вычислительного устройства для оптимального выбора параметров при распараллеливании. Под этим подразумевается количество доступных вычислительных ядер на многопроцессорном вычислительном устройстве или графической карте. В данной работе не применялись специализированные математические библиотеки или другие возможности и распараллеливание осуществлялось разбиением двухмерных или одномерных циклов «for» при обходе 1D границ или 2D точек слоев. Отметим, что программа расчетов реализована с применением неблокированных функций обмена данными по MPI. Такая организация позволят выполнить пересылки за то время, что происходит расчет внутренней двухмерной части каждого из слоев. Таким образом вычислительные процедуры представлены несколькими типами: одни - для расчета точек внутри слоя на основе 2D процедур, другие - для расчета точек на смежных границах соседних слоев. Поэтому и распараллеливание применятся либо по двум координатам, либо только по одной координате. В случае использования в расчетах GPU не стоит забывать о необходимости копирования данных между CPU и GPU, которые реализуются с использованием специализированных буферов данных, инициализированных как на управляющем
CPU, так и на вычислительном устройстве GPU. Более подробно о таком способе организации вычислений можно ознакомиться в работе [30], где аналогичная технология применялась для решения динамической задачи теории упругости для трехмерного варианта на основе конечно-разностного метода второго порядка аппроксимации. В отличие от параллельного алгоритма для разностной схемы второго порядка в реализации для схемы четвертого порядка необходимо передавать по два слоя точек вдоль оси Oz «сверху» и «снизу» каждого слоя, исключая самый верхний и самый нижний в топологии, рис. 3.
Рис. 3. Схематическое представление организации параллельных вычислений для предложенного алгоритма
Расчеты и результаты
Разработанная технология математического моделирования и созданные программы предназначены для возможности расчетов сейсмического поля на большие расстояния. Это связано с исследованием строения и определением значений упругих параметров для Байкальской рифтовой зоны. Такой объект характеризуется большой протяженностью. Для подготовки сеточной модели использовались данные, которые содержат информацию о значениях скоростей упругих волн и плотности на выделенных горизонтах. Было проведено восстановление распределения значений упругих партеров на расчетной сетке на основе разработанной программы для создания моделей упругих сред. В такой программе применяется алгоритм сплайн интерполяции. Результатом работы стало восстановление трех массивов данных для проведения расчетов и определения значений X и ц. Область моделирования составила 400 км по оси Ox и 73 км по оси Oz, рис. 4. Для восстановления использовано более 10 опорных точек (скважинных данных).
Для расчета сейсмического поля в программе использовался источник типа «вертикальная сила» с частотой 7 Гц. Координаты расположения источника ~15 км по оси Ox. Расчётная сетка включала 56743 х 10356 точек. Размер одного выходного файла с результатами в виде снимка волнового поля составил более 2 Гб. Расчет проведен на кластере НКС-30Т Центра коллективного пользования Сибирского суперкомпьютерного центра (http://www.sscc.icmmg.nsc.ru). Для расчетов применена программа на основе работы с CPU и MPI. В расчетах
использовано 12 вычислительных узлов по 12 ядер. На каждое ядро выделено по одному MPI процессу.
Рис. 4. Геометрическое строение модели Байкальской рифтовой зоны на основе скважинных данных. Рисунок создан по распределению значений скорости S волн на расчетной сетке
Рис. 5. Снимки вертикальной компоненты сейсмического поля рифтовой зоны для времен расчета: 15, 20, 25, 30, 40 секунд (сверху-вниз)
Заключение
Предложен алгоритм решения динамической задачи теории упругости для моделирования полного поля сейсмических волн от точечного источника. Предложена параллельная реализация конечно-разностного метода для численных экспериментов. По результатам работы разработана технология математического моделирования и программы для проведения расчетов на высокопроизводительных вычислительных системах с вычислительными устройствами CPU и GPU. С использованием разработанных программ подготовлена геофизическая модель строения Байкальской рифтовой зоны для расстояний на сотни километров. Для восстановления структуры и распределения упругих параметров на сеточной модели использовались скважинный данные, содержащие информацию о значениях на выделенных горизонтах. Для подготовленной модели выполнено математическое моделирование процессов распространения сейсмических волн для профиля на 400 км. Для данного объекта исследований получены синтетические данные в виде сейсмограмм и снимков волнового поля. С использованием полученных результатов возможно провести исследовательскую работу по сопоставлению синтетических и экспериментальных данных, что позволит уточнить и скорректировать геофизическую модель.
Работа выполнена при поддержке РФФИ (гранты № 16-07-01052, 17-07-00872).
REFERENCES
1. Maeda T., Takemura S. and Furumura T. (2017). OpenSWPC: an open-source integrated parallel simulation code for modeling seismic wave propagation in 3D heterogeneous viscoelastic media. Earth, Planets and Space 69:102 DOI 10.1186/s40623-017-0687-2, https://link.springer.com/content/pdf/10.1186/s40623-017-0687-2.pdf.
2. Komatitsch, D., et al. (2010). High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster. J. Comput. Phys. 229(20), 7692-7714.
3. Aochi H., Ulrich T., Ducellier A., Dupros F., Michea D. (2012) Finite difference simulations of seismic wave propagation for understanding earthquake physics and predicting ground motions: Advances and challenges. Conference on Computational Physics. https://hal-brgm.archives-ouvertes.fr/hal-00726430/document.
4. Mordvinova V.V., Artemiev A.A. (2010). The three-dimensional shear velocity structure of lithosphere in the southern Baikal rift system and its surroundings. Russian Geology and Geophysics, Vol. 51, No. 6, pp. 694-707.
5. O'Reillya O., Lundquist T., M.Dunhama E., Nordstromb J. (2017). Energy stable and high-order-accurate finite difference methods on staggered grids Journal of Computational Physics 346 572-589.
6. Solano P. C.A., Donno D., Chauris H. (2016). Finite-difference strategy for elastic wave modelling on curved staggered grids, Comput. Geosci. 20(1) 245-264, http://dx.doi.org/10.1007/s10596-016-9561-8.
7. O'Reilly O., Nordstrom J., Kozdon J.E., Dunham E.M. (2015). Simulation of earthquake rupture dynamics in complex geometries using coupled finite difference and finite volume methods, J. Commun. Phys. 17(2) 337-370, http://dx.doi.org/10.4208/cicp.111013.120914a.
8. Ren Z., Liu Y. (2015) Acoustic and elastic modeling by optimal time-space-domain stag-gered-grid finite-difference schemes, Geophysics, 80(1), T17-T40.
9. Yong P., Huang J., Li Z., Liao W., Qu L., Li Q., Liu P. (2017). Optimized Equivalent Staggered-grid FD Method for Elastic Wave Modeling Based on Plane Wave Solutions, Geophysical Journal International, Volume 208, Issue 2, pp. 1157-1172.
10. Heinecke A., Breuer A., Bader M., Dubey P. (2016). High Order Seismic Simulations on the Intel Xeon Phi Processor (Knights Landing). High Performance Computing, pp 343-362.
11. Tobin J., Breuer A., Heinecke A., Yount C., Cui Y. (2017) Accelerating Seismic Simulations Using the Intel Xeon Phi Knights Landing Processor. In: Kunkel J., Yokota R., Balaji P., Keyes D. (eds) High Performance Computing. ISC 2017. Lecture Notes in Computer Science, vol 10266. Springer. pp 139-157.
12. Li Y., Metivier L., Brossier R., Han B. and Virieux J. (2015). 2D and 3D frequency-domain elastic wave modeling in complex media with a parallel iterative solver Geophysics, Vol. 80, No. 3; P. T101-T118.
13. Chan J., Wang J., Modave A., Remacle J.-F., and Warburton T. (2016). GPU-accelerated discontinuous Galerkin methods on hybrid meshes. Journal of Computational Physics Volume 318, Pages 142-168.
14. Castro M., Francesquini E., Dupros F., Aochi H., Navaux P., et al. (2016). Seismic Wave Propagation Simulations on Low-power and Performance-centric Manycores. Parallel Computing, Elsevier, Volume 54, Pages 108-120.
15. Qawasmeh A., Chapman B., Hugues M., Calandra H. (2015). GPU Technology Applied to Reverse Time Migration and Seismic Modeling via OpenACC. Proceeding PMAM '15 Proceedings of the Sixth International Workshop on Programming Models and Applications for Multicores and Manycores Pages 75-85.
16. Iturraran-Viveros U. and Molero-Armenta M. (2015). GPU computing with OpenCL to model 2D elastic wave propagation: exploring memory usage Comput. Sci. Disc. 8.
17. Hamilton B., Webb C. J., Gray A., Bilbao S. (2015). Large Stencil Operations for GPU-based 3-D Acoustics Simulations Proc. of the 18th Int. Conference on Digital Audio Effects (DAFx-15), https://www.ntnu.edu/documents/1001201110/1266017954/DAFx-15_submission_46.pdf.
18. Darmawan J. B. B. and Mungkasi S. (2017). Performance of parallel computation using CUDA for solving the one-dimensional elasticity equations J. Phys.: Conf. Ser. 801 012080, doi:10.1088/1742-6596/801/1/012080, http://iopscience.iop.org/article/10.1088/1742-6596/801/1/ 012080/pdf.
19. Assi H., Cobbold R. S. (2016). A second-order, perfectly matched layer formulation to model 3D transient wave propagation in anisotropic elastic media. The Journal of the Acoustical Society of America 140, 3261 https://doi.org/10.1121/L4970329.
20. Lee and C. Shin. (2015). Time-domain formulation of a perfectly matched layer for the second-order elastic wave equation with VTI media. J. Seism. Explor, vol. 40, pp. 231-257.
21. Xie Z., Matzen R., Cristini P., Komatitsch D., and Martin R. (2016). A perfectly matched layer for uidsolid problems: Application to ocean-acoustics simulations with solid ocean bottoms. The Journal of the Acoustical Society of America, 140(1):165-175.
22. Berenger J.-P. (1994). A perfectly matched layer for the absorption of electromagnetic waves. J. Comp. Phys. 114, 185-200.
23. Collino F. and C. Tsogka. (1998). Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media. INRIA, Rapportde recherche, n. 3471.
24. Chew W. C. and Q. H. Liu. (1996). Perfectly matched layers for elastodynamics: a new absorbing boundary condition. J. Comp. Acoustics 4, 341-359.
25. Komatitsch, D. and J. Tromp. (2003). A Perfectly Matched Layer (PML) absorbing condition for the second-order elastic wave equation. Geophys. J. Int. 154, 146-153.
26. Komatitsch D., Martin R. (2007). An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation. Geophysics, vol. 72, pp. 155-167.
27. Levander A. (1988). Fourth-order finite difference P-SV seismograms, Geophysics 53, 1425-1436.
28. Moczo P. (1998). Introduction to Modeling Seismic Wave Propagation by the Finite-Difference Method. Lecture Notes, Kyoto University, ftp://www.nuquake.eu/pub/Papers/ Moczo_LN_Kyoto_1998.pdf.
29. CUDA C Programming Guide. (2017). http://hpc.pku.edu.cn/docs/ 20170829223400085039.pdf.
30. Karavaev D.A., Glinsky B.M., Kovalevsky V.V. (2015). A Technology of 3D Elastic Wave Propagation Simulation Using Hybrid Supercomputers. CEUR Workshop Proceedings 1st Russian Conference on Supercomputing Days, Vol. 1482, pp. 26-33.
© Д. А. Караваев, 2018
