CC BY
24Математические заметки СВФУ Июль—сентябрь, 2015. Том 22, № 3
УДК 519.688
ЭЭ-МОДЕЛИРОВАНИЕ СЕЙСМИЧЕСКИХ ВОЛНОВЫХ ПОЛЕЙ В СРЕДАХ, ХАРАКТЕРНЫХ ДЛЯ ВУЛКАНИЧЕСКИХ СТРУКТУР
Б. М. Глинский, В. Н. Мартынов, А. Ф. Сапетина
Аннотация. Необходимость предсказания катастрофических явлений, которые могут быть вызваны готовящейся вспышкой вулканической деятельности, является актуальной задачей. Для решения этой задачи необходимо проведение комплексных и объективных исследований процессов, происходящих на поверхности и внутри вулканической структуры, в том числе математическое моделирование с целью создания системы вибросейсмического мониторинга. Разработаны параллельные 2D- и SD-алгоритмы и созданы программы для моделирования распространения упругих волн в сложно построенной SD-упругой среде ^D-модель есть сечение исходной SD-модели различными плоскостями и под разными углами). Используется конечно-разностный метод, явные разностные схемы на сдвинутых сетках и метод поглощающих границ CFS-PML. Предложенный численный метод и его параллельная программная реализация эффективно используют архитектуру современного суперкомпьютера на основе графических ускорителей. Проведена серия SD-расчетов с целью изучения структуры волнового поля, обусловленного геометрией внутренних границ, для уточнения кинематических и динамических характеристик волнового поля с целью возможного создания системы вибросейсмического мониторинга стратовулкана Эльбрус.
Ключевые слова: мониторинг, SD-моделирование, упругие волны, разностные схемы, гибридный кластер, GPU.
B. M. Glinskiy, V. N. Martynov and A. F. Sapetina. 3D-Modeling of Seismic Wave Fields in a Medium Specific to Volcanic Structures. Abstract: The need for prediction of catastrophic events that may be caused by the impending outbreak of volcanic activity, is an actual problem. It is necessary for solve this problem to conduct comprehensive and objective investigation of the processes occurring on the surface and inside the volcanic structure, including mathematical modeling to create a system of vibroseis monitoring. 2D and 3D-algorithms are developed and software for simulation of elastic wave propagation in a complicated medium (2D model is separation of original 3D-model using various angles and planes) is created. Explicit finite-difference schemes for the shifted grids and CFS-PML method of absorbing boundaries are used. Proposed numerical method and its parallel realization efficiently uses architecture of modern supercomputer equipped with graphical accelerator. Serial of 3D-estimation is carried out for the structure researching of wave field due to the geometry of the internal boundaries and refinement of its kinematic and dynamic characteristics in order to create the system of vibroseis monitoring of stratovolcano Elbrus. Keywords: monitoring, 3D-simulation, elastic waves, finite difference schemes, hybrid cluster, GPU.
© 2015 Глинский Б. М., Мартынов В. Н., Сапетина А. Ф.
1. Введение
Современный процесс разработки программного инструментария для проведения численного моделирования становится все более зависимым от архитектуры суперкомпьютера, подчиняя себе не только выбор алгоритмов для решения поставленных задач под вычислительную архитектуру, но и выбор самой постановки. Разные алгоритмы быстрее работают на разных типах вычислительных архитектур и практически всегда требуется некая адаптация выбранного алгоритма под имеющуюся архитектуру, которая может дать большой прирост производительности при решении задачи, что соответственно уменьшает время моделирования. Данная работа опирается именно на такой подход к разработке инструментария для численного моделирования сейсмических волновых полей в 3Б-средах, характерных для вулканических структур.
Вулканы, в том числе и спящие, являются потенциальной угрозой внезапного возникновения мощных катастрофических явлений. Возможность предсказать заранее приближающееся извержение позволила бы сохранить жизни людей и их имущество. Для этого необходимо проводить мониторинг состояния вулкана, отслеживать изменения параметров внутренней среды, их вариации, и делать по ним прогноз извержений. Последнее требует предварительного тщательного комплексного исследования процессов, проходящих внутри вулканов и на их поверхности.
Одним из инструментов для проведения такого исследования является активный вибросейсмический мониторинг. Получаемая в его ходе информация сложна для интерпретации и требует проведения предварительного и сопутствующего математического моделирования процессов, происходящих внутри исследуемого геофизического объекта с учетом особенностей его строения. Обычно рельеф исследуемого объекта достаточно сложен и не позволяет расположить на нем площадную систему наблюдения для решения обратной задачи геофизики, поэтому приходится решать прямую задачу, варьируя параметры моделируемой среды так, чтобы результаты численного и натурного экспериментов совпадали.
Кроме того, такой подход позволяет рассмотреть возможные типы строения магматических вулканов и их различные состояния, чтобы выделить основные эффекты, возникающие в данных, фиксируемых системой наблюдения на поверхности вулкана, для облегчения их интерпретации при непосредственном мониторинге.
Таким образом, мы сталкиваемся с проблемой проведения крупномасштабного численного моделирования процессов, происходящих в вулканических структурах различного строения при вибросейсмическом мониторинге. Решение такой задачи требует использования новейших суперкомпьютерных технологий.
2. Постановка задачи и метод решения
Численное моделирование распространения сейсмических волн в упругих неоднородных средах проводится на основе решения полной системы уравнений
теории упругости с соответствующими начальными и граничными условиями, записанной в терминах вектора скоростей смещений и = (и, V, Ш)т и тензора напряжений а = (ахх ,ауу , аху, ахг, аух)т.
В качестве области моделирования рассматривается изотропная ЗБ-неод-нородная сложно построенная упругая среда, представляющая собой параллелепипед, одна из граней которого является свободной поверхностью (плоскость 2 = 0).
Основные уравнения в векторной форме могут быть представлены в следующем виде:
ди , -а, ч да ,
Р^Т = [А]ст + Р(г,х,у,г), 1- = \В\и,
В
г — дх 0 0 # ду д дх 0
А = 0 _д_ ду и дх 0 д дх ,
0 0 £ 0 д_ дх д_ ду -
Г (А + -2м) _д_ дх ду дх
(А+ 2//)^ А# дх
ду (А+ 2
Рдх 0
0
_ 0 »Ш <
где Ь — время, р(х,у,г) — плотность, Х(х,у,г), ^(х,у,г) Начальные условия имеют вид
(1)
параметры Ламе.
ахх |^=0 = 0, аух |^=0 = 0, аху ^=0 = 0, ахх^=0 = ауу ^=0 = ахх ^=0 = 0,
(2)
и |,=0 = 0, V |4=0 = 0, Ш |4=0 = 0,
а граничные условия на свободной поверхности — вид
аХХ |2 = 0 = 0 ауг |2 = 0 = 0 ахх |2 = 0 = 0- (3)
Для численного решения поставленной задачи (1)-(3) применяется известная конечно-разностная схема Верье [1-3], хорошо себя зарекомендовавшая. Расчет ее сеточных коэффициентов проводится на основе интегральных законов сохранения. Схема имеет второй порядок аппроксимации по времени и пространству [1], в данной работе рассматриваются только равномерные сетки. Для примера представим несколько конечно-разностных уравнений используемой схемы:
„.п+1 „,п п+± п+±
„ I/-, и 1-, — и 1-7 /Т 2 _ /Т 2
Рг,],к т _ ихх1,],к ихх1-1^,к
2 Т ~ Дж
а т-.^т а 1 . 1, а 1-711 а 1-7 1
хуг-±,з+±,к хуг-±,з-±,к хгг-±,з,к+± хгг-±,з,к-± п
Ау Аг ^ !хг^
п+4 п—4
(7ХХЪ_ \ 2 /с-1 ~~ ^ А;--
= /^-им -к;-+-д^-
где
! 1 ! 1 1 1 1 \\ -= (т(-+-+-+-))
4 Мг —1 — 1 Мг——1
Границы области моделирования вызывают ложные отражения внутри нее. Для поглощения этих отражений используется вспомогательный метод СЕ8-РМЬ [4-6], имеющий ряд преимуществ перед классическим методом РМЬ. Он дает более «качественную» картину волнового поля для данной задачи, более прост в реализации и экономичен с вычислительной точки зрения. Для применения метода каждая из границ моделируемой области — параллелепипеда — за исключением свободной поверхности, располагающейся на верхней грани, окружается поглощающим слоем. Во внутренности волновое поле рассчитывается по первоначальным конечно-разностным уравнениям, а при попадании волны в зону поглощения происходит расчет по новым формулам с демпфирующими параметрами, описывающими подход к созданию поглощающих границ. Выбор значений демпфирующих параметров для численных расчетов в соответствующих поглощающих слоях проводится на основе результатов работы [5].
3. Программный комплекс для проведения численных экспериментов
В данной работе предлагается подход для численного моделирования распространения сейсмических волн в средах, характерных для вулканических структур, в рамках которого требуется разработать комплекс алгоритмов и программ, позволяющих проводить конструирование сеточной модели среды сложной геометрии и необходимые расчеты. При этом предлагается в зависимости от текущей цели и имеющихся ресурсов проводить или ЗБ-моделирование во всей интересующей области, или 2Б-моделирование в сечениях изначальной области, содержащих наиболее интересные характерные особенности волнового поля, или сочетать эти подходы.
Таким образом, предлагаемый программный комплекс должен состоять из следующих частей:
— программа для конструирования сеточных моделей сложно построенных сред с включениями, характерными для магматических вулканов;
— программа численного моделирования распространения упругих волн в ЗБ-неоднородных упругих средах с криволинейной свободной поверхностью;
— программа численного моделирования распространения упругих волн в ЗБ-неоднородных упругих средах с прямолинейной свободной поверхностью;
— программа численного моделирования распространения упругих волн в 2Б-неоднородных упругих средах с криволинейной свободной поверхностью для заданного сечения рассматриваемой ЗБ-модели;
— программа численного моделирования распространения упругих волн в 2Б-неоднородных упругих средах с прямолинейной свободной поверхностью для заданного сечения рассматриваемой ЗБ-модели.
На данный момент реализованы части предлагаемого комплекса, позволяющие работать со случаем прямолинейной свободной поверхности. Разработка программ велась с учетом особенностей архитектуры гибридного кластера НКС-ЗОТ+GPU ССКЦИВМиМГ СО РАН (http://www2.sscc.ru), включающего в себя 40 вычислительных узлов, оснащенных графическими картами NVIDIA Tesla M2090 на архитектуре Fermi. Пиковая производительность 85 Тфлопс.
Для эффективного использования гибридной архитектуры необходима адаптация и оптимизация используемых в ходе моделирования алгоритмов, основанная на знании архитектуры кластера и его составляющих и необходимых программных средств. Подробное описание разработанных программных средств и параллельной реализации, адаптации и оптимизации алгоритмов приведено в [7,8]. Повторим здесь некоторые основные моменты.
В разработанный комплекс входит построитель моделей среды, который позволяет конструировать на сеточном уровне сложные модели упругих сред на основе идеи Z порядка, описанной в [2,3]. При параллельном исполнении данные сразу конструируются на тех вычислительных узлах, на которых впоследствии будет проводиться расчет.
Также создана параллельная программа для численного моделирования распространения волн в трехмерных неоднородных упругих средах с прямолинейной свободной поверхностью, реализующая указанные ранее методы (схема Верье и CFS-PML), адаптированные под архитектуру гибридного кластера.
Для распараллеливания данной задачи используется декомпозиция области на слои вдоль направления одной из координатных осей. Каждый слой рассчитывается на отдельном узле, где, в свою очередь, он разбивается еще на подслои вдоль другой координатной оси по числу графических ускорителей на узле. При такой реализации каждая графическая карта рассчитывает свою сеточную область внутри подслоя на каждом временном шаге независимо от других, за исключением точек, находящихся на границе между двумя соседними областями. Эти точки являются общими для каждой из областей, и для продолжения счета необходимо производить обмен информацией об искомых величинах между графическими картами на узле и между соседними узлами в направлении разных координатных осей. Обмены производятся при помощи технологии MPI (Message Passing Interface), работа с графическими ускорителями осуществляется с помощью технологии CUDA (Compute Unified Device Architecture). При этом на каждой графической карте параллельная часть кода выполняется как большое количество нитей. Такой гибридный подход обеспечивает высокую степень параллелизма.
Также в разработанный комплекс входит программа для численного моделирования распространения сейсмических волн в двумерной упругой среде, которая разработана на основе программы для трехмерного случая. Для решения 2Б-задачи используется та же разностная схема на сдвинутых сетках и вспомогательный метод CFS-PML в их двумерном варианте.
Для 2Б-расчета требуется намного меньше вычислительных ресурсов, в том числе и памяти для хранения всех необходимых данных, поэтому данная реализация использует всего один вычислительный узел гибридного кластера с тремя графическими ускорителями. В итоге разработанная программа позволяет проводить полномасштабные расчеты за приемлемое время, задействовав при этом намного меньшее количество ресурсов по сравнению с 3Б-вариантом, который потребует для расчета решения практически все ресурсы гибридного кластера и намного больше времени для вычислений. Для примера, построение модели и реальный расчет волнового поля на узле с тремя графическими ускорителями для соответствующих сеток по времени и пространству (6000 х 9000 узлов по пространству и 25000 шагов по времени) занимает всего 12 минут.
4. Исследование масштабируемости разработанного ПО
Для анализа производительности разработанного программного обеспечения исследована его сильная и слабая масштабируемость, проведена имитация исполнения программы на большом количестве ядер и сравнение времени работы программы на гибридном кластере с аналогичным временем для программы, исполнявшейся на кластере с классической MPP-архитектурой.
Под сильной масштабируемостью будем понимать уменьшение времени счета одного шага одной и той же задачи при использовании большего числа графических ядер в рамках одной графической карты. Ее исследование позволяет понять, насколько эффективно применяемые алгоритмы используют архитектуру самой графической карты. Под слабой масштабируемостью будем понимать сохранение времени счета одного шага одного и того же объема задачи при одновременном увеличении количества графических карт.
Результаты исследований для программы 3Б-моделирования приведены на рис. 1 и 2 в виде графиков. Из графика на рис. 1 видно, что задача хорошо ложится на архитектуру графической карты (получено ускорение около 40 раз при использовании всех ядер GPU по сравнению с одним ядром GPU).
Из графика на рис. 2 видно, что достигнута эффективность около 92% при увеличении количества графических карт до 30 штук. На основе собранных данных о времени вычисления искомых компонент тензора напряжения и вектора скоростей смещения и времени обменов между вычислительными узлами и графическими ускорителями для реализованного программного обеспечения была проведена имитация выполнения алгоритма численного моделирования распространения сейсмических волн в упругой среде с помощью распределенной агентно-ориентированной системы имитационного моделирования AGNES, разработанной в ИВМиМГ.
Пакет AGNES (AGent NEtwork Simulator) базируется на Java Agent Development Framework (JADE). JADE — это мощный инструмент для создания мультиагентных систем на JAVA, и он состоит из 3-х частей: среда исполнения агентов; библиотека базовых классов, необходимых для разработки агент-ной системы; набор утилит, позволяющих наблюдать и администрировать МАС (мультиагентная система). Для моделирования больших вычислений важно,
Рис. 1. График исследования сильной масштабируемости.
Рис. 2. График исследования слабой масштабируемости.
что JADE — это FIPA-совместимая распределенная агентная платформа, которая может использовать один или несколько компьютеров (узлов сети), на каждом из которых должна работать только одна виртуальная JAVA машина [9]. Данная система позволяет отображать вычислительный алгоритм на гипотетический суперкомпьютер, исследовать его поведение, провести корректировку вычислительной схемы. Примеры применения данного подхода изложены в [10].
Результаты имитации представлены на рис. 3. Начало имитационного процесса сопоставлено с данными реального исследования слабой масштабируемости. Имитация показывает, что при увеличении количества вычислительных ядер до миллиона эффективность работы предлагаемой программы будет равна около 75%. Это говорит о том, что разработанное ПО обладает хорошей масштабируемостью, что позволит использовать его на больших вычислительных комплексах для проведения крупномасштабных исследований.
Проведено сравнение времени, затраченного на численное моделирование 3D-модели с использованием для вычислений узлов с GPU, с аналогичным временем для расчета на классическом кластере с CPU. Для этого используется время полномасштабного расчета, указанное в [2] для сетки размером 1677 х 1059 х 971 по пространству и 10313 шагов по времени. Расчет проведен
Рис. 3. Зависимость относительного ускорения SL(M)
от общего числа моделируемых GPU ядер M (горизонтальная ось — в логарифмическом масштабе)
на 20 вычислительных блэйд-серверах hp ProLiant BL2x220c G5, находящихся в составе НКС-30Т ССКЦ ИВМиМГ СО РАН. Соответствующее время составило 31 час 15 мин. 17 с (задействовано 160 CPU ядер). Аналогичное время расчета с теми же размерами расчетной сетки на 15 узлах гибридного кластера с помощью разработанного программного обеспечения составило 2 час 56 мин. (задействовано 15360 GPU ядер); полученное ускорение составило 10,66 раза.
Таким образом, созданное программное обеспечение, адаптированное под современную гибридную вычислительную архитектуру с графическими ускорителями, использующимися в различных вычислительных центрах России и мира, обладает высокой производительностью.
5. Численные эксперименты
На основании данных, приведенных в [11-15], а также в работах других авторов, построена геофизическая модель стратовулкана Эльбрус. Вулканическая постройка лежит на гранитном блоке +I (рис. 4); эффузивные породы слагают вулканический конус +II; ниже нулевой отметки можно выделить 8 слоев (табл. 1). Зададим верхний магматический очаг в виде эллипсоида с горизонтальными и вертикальными осями 9 и 6 км (р = 2,1г/см3, Vp = 2, 2 км/с); диаметр бывшего канала 130 м; зададим материнский магматический очаг: эллипсоид с горизонтальными и вертикальными осями 24 и 13 км (р =1, 8 г/см3, Vp = 1, 9 км/с), диаметр предполагаемого подпитывающего канала 250 м. Средний канал — цилиндр диаметром 160 м. Итак, в качестве приближенной модели вулкана Эльбрус можно принять многослойную среду с включениями в виде эллипсов цилиндров с параметрами, указанными в табл. 1. Подробное описание геофизической модели можно найти в [8].
Построенная модель была взята за основу для проведения последующих вычислительных экспериментов, результаты которых приведены далее для иллюстрации работоспособности разработанного комплекса программ, а также для
Рис. 4. Геофизическая модель вулкана Эльбрус и схема вибросейсмического мониторинга.
Таблица 1. Параметры среды для геофизической модели вулкана Эльбрус
Vp, км/с Vs, км/с р, г/см3
Слой +11 2,85 1,65 2,4
Слой +1 3,1 1,79 2,66
Слой I 3,2 1,82 2,7
Слой II 5,9 3,42 2,85
Слой III 6,22 3,59 2,62
Слой IV 5,82 3,37 2,7
Слой V 5,97 3,45 2,75
Слой VI 6,43 3,72 2,78
Слой VII 6,95 4,03 2,81
Слой VIII 8,1 4,68 2,85
того чтобы продемонстрировать различные возможности применения предложенной технологии суперкомпьютерного моделирования распространения сей-
смических волн в средах характерных для вулканических структур.
Каждый расчет проводился на 11 узлах кластера НКС-30Т+СРи, оснащенных графическими картами. Для моделирования процесса на 12000 шагов по времени тратится в среднем чуть больше полутора часов.
Результаты 3Б-численного моделирования содержат большой объем информации. В связи с этим для представления и анализа результатов численного моделирования мы используем теоретические сейсмограммы и мгновенные снимки различных сечений 3Б-волнового поля плоскостями, проходящими через линии возможной системы наблюдений.
Система возбуждения для всех расчетов состоит из точечного источника типа «центр давления» с частотой 8 Гц, располагающегося вблизи свободной поверхности в левой части расчетной области в одной плоскости с осью симметрии магматических каналов и камер.
Для первого расчета в качестве среды взят фрагмент исходной приближенной модели стратовулкана Эльбрус, включающий в себя только верхнюю магматическую камеру видоизмененной формы и прилегающие к ней каналы, расположенные в пятислойной среде. Параметры среды взяты из табл. 1. Целью расчета было проиллюстрировать возможности построителя сеточных моделей среды, с помощью которого можно сконструировать среду с включениями достаточно сложной формы (например, пересечение нескольких объектов различной природы под различными углами). В данном случае магматическая камера имеет форму двух пересекающихся эллипсоидов.
Результаты расчета приведены на рис. 5 в виде мгновенных снимков волнового поля в плоскости XZ, проходящей через точечный источник и ось симметрии верхнего канала. Визуализация снимков проведена с помощью программы Лвр1в, разработанной в ОАО «Сибнефтегеофизика». На снимках отмечены границы между различающимися по параметрам участками моделируемой среды.
Из рис. 5 видно, что волновое поле имеет сложную структуру и существенно зависит от геометрии, размеров и свойств включений. Для проведения геофизической интерпретации теоретических сейсмограмм, полученных также в ходе расчета такой сложно построенной среды, требуется проведение большой серии вычислительных экспериментов для 3Б-модели среды, а также ее характерных 2Б-сечений, которые могу быть использованы с этой целью.
В [8] приведено сравнение результатов 3Б- и 2Б-расчетов для фрагмента предлагаемой приближенной модели стратовулкана Эльбрус. Представленные снимки являются иллюстрацией нашего предположения о том, что 2Б-моделирование может быть использовано для возможного построения 3Б-моде-лей среды. Расчет может производиться одновременно для 3Б-модели среды и нескольких 2Б-сечений для получения необходимой информации об особенностях строения волнового поля для выбранной геофизической модели, системы наблюдения и положения источника. Такой подход значительно ускоряет процесс исследования особенностей волнового поля при моделировании процесса вибросейсмического зондирования вулканической постройки и позволяет провести меньшее количество расчетов для трехмерной задачи, требующих значительно большего числа ресурсов.
Т = 1,15 с
Т = 1,84 с
Т = 2, 53 с
Т = 3, 22 с
Рис. 5. Мгновенные снимки сечения 3Б-волнового поля в различные моменты времени для пятислойной среды с включением из двух пересекающихся эллипсоидов (компонента и, плоскость XZ).
())(*))$%($(')#$%%(+*$((#))+%'(%%)+#(+*%(**+)&'(*%$%*++ )**&%)%(&((')&*%(&(&&)*&($*+)%'&#)$#&&)%)#((
01 X = 2,07 е
02 X = 2,76 е
Рис. 6. Мгновенные снимки сечения ЗБ-волнового поля в различные моменты времени (компонента и, плоскость XZ): А — среда без верхнего канала (равновесное состояние вулкана), В — среда, содержащая верхний канал, заполненный магмой (момент извержения), С — разница между сейсмограммами А и В.
Разработанный комплекс программ создавался с целью проведения исследований процессов, проходящих внутри магматических вулканов и на их поверхности, для осуществления последующего мониторинга извержений. Первым шагом на пути к этой цели могут быть следующие расчеты.
В качестве моделируемой среды для них также взят фрагмент предложенной приближенной модели стратовулкана Эльбрус, включающий в себя только верхнюю магматическую камеру в форме эллипсоида и прилегающий к ней подпитывающий канал, расположенные в пятислойной среде. Разница между этими средами заключается в том, что в одном случае верхний канал заполнен магмой (момент извержения), а во втором случае он сливается с окружающими слоями (равновесное состояние). Параметры среды взяты из описания выше.
Результаты расчетов представлены на рис. 6 в виде мгновенных снимков волнового поля в плоскости XZ, проходящей через точечный источник и ось симметрии верхнего канала. Также добавлены разницы между волновым полем соответствующих снимков. На них наглядно видны эффекты, возникающие в связи с присутствием верхнего канала.
Отметим важность проведения интерпретации теоретических сейсмограмм, получаемых в ходе численного эксперимента, так как в ходе натурных экспериментов исследователи получают данные именно с поверхности, соответственно сопоставление теоретических и натурных данных может быть проведено именно на этом уровне. При этом снимки волнового поля, полученные в ходе расчета,
ABC
Рис. 7. Теоретические сейсмограммы для компоненты U: A — среда без верхнего канала (равновесное состояние вулкана),
B — среда, содержащая верхний канал, заполненный магмой (момент извержения), С — разница между сейсмограммами A и B.
могут служить материалом для объяснения эффектов, возникающих на сейсмограммах.
На рис. 7 приведены теоретические сейсмограммы для проведенных расчетов, а также их разница, на которой видны возникающие различия в общей картине.
Полученные данные сложны для интерпретации. Тщательное изучение теоретических сейсмограмм на основе полученных снимков волнового поля и дополнительных расчетов может выявить признаки, по которым будет осуществляться мониторинг извержений.
Заключение
В работе описана технология решения задачи распространения упругих волн в сложно построенных средах, учитывающая различные по размерности численные методы исследования неоднородной среды и ориентированная на гибридную архитектуру суперкомпьютера, оснащенного графическими ускорителями.
Разработан комплекс параллельных программ с учетом выбранной архитектуры, реализующий распространение упругих волн в 2D- и 3Б-моделях в сложно построенных средах, характерных для вулканических структур. Проведено исследование времени работы программного обеспечения, в том числе имитация выполнения программы на большом числе ядер, которое показывает, что созданное ПО обладает высокой производительностью и может эффективно
использовать ресурсы больших вычислительных комплексов.
Различные численные эксперименты для различных вариаций предложенной упрощенной модели вулкана Эльбрус показывают работоспособность созданного комплекса программ и различные возможности его применения для последующих исследований. Во всех экспериментах рассчитанное волновое поле от точечного источника имеет сложную картину и существенно зависит от геометрии, размеров и свойств включений. Интерпретация данных, полученных в ходе численного моделирования для сложно построенных сред с наличием различных включений, обусловленных приближением к реальной модели среды, — непростая задача. Тем не менее представленные эксперименты показывают, что численное моделирование может дать значительную информацию для организации проведения натурных экспериментов и интерпретации результатов экспериментальных наблюдений, полученных в ходе вибросейсмического мониторинга. Особенности, которые можно выявить из полученных результатов, могут быть использованы при вибросейсмическом мониторинге вулканических структур. Проведение серий подобных расчетов позволяет выбрать модель вулканической структуры, адекватную результатам натурных наблюдений по кинематическим и динамическим характеристикам.
В дальнейшем предполагается развитие предложенной суперкомпьютерной технологии численных расчетов и последующее проведение численных экспериментов для моделей с криволинейной свободной поверхностью.
ЛИТЕРАТУРА
1. Bihn M., Weiland T. A. Stable discretization scheme for the simulation of elastic waves // Proc. 15th IMACS World Congress Scientific Computation, Modelling and Applied Mathematics (IMACS 1997). 1997. V. 2. P. 75-80.
В 1 надо указать место издания и издательство.
2. Глинский Б. МКараваев Д. А., Ковалевский В. В., Мартынов В. Н. Численное моделирование и экспериментальные исследования грязевого вулкана «Гора Карабетова» вибросейсмическими методами // Вычисл. методы и программирование. 2010. Т. 11. С. 95-104.
3. Караваев Д. А. Параллельная реализация метода численного моделирования волновых полей в трехмерных моделях неоднородных сред // Вестн. Нижегород. ун-та им. Н. И. Лобачевского. 2009. Т. 6, № 1. С. 203-209.
4. Drossaert F., Giannopoulos A. Complex frequency shifted convolution PML for FDTD modeling of elastic waves // Wave Motion. 2007. V. 44. P. 593-604.
5. Komatitsch D., Martin R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation // Geophysics. 2008. V. 73, N 4. P. T51-T61.
6. Hastings F. D., Schneider J. B., Broschat S. L. Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation // J. Acoust. Soc. Am. 1996. V. 100, N 5. P. 3061-3069.
7. Сапетина А. Ф. Численное моделирование распространения сейсмических волн в сложно построенных средах на гибридном кластере // Пробл. прочности и пластичности. 2014. Т. 76, № 4. С. 288-296.
8. Глинский Б. М., Мартынов В. Н., Сапетина А. Ф. Технология суперкомпьютерного 3D-моделирования сейсмических волновых полей в сложно построенных средах // Вестн. Южно-Урал. гос. ун-та. Сер. Вычисл. математика и информатика. 2015. Т. 4, № 4. С. 101-116.
9. Подкорытов Д. И. Агентно-ориентированная среда моделирования сетевых систем AGNES // Ползунов. вестн. 2014. № 2/1. С. 93-106.
10. Глинский Б. М., Куликов И. М., Снытников А. В., Черных И. Г., ВинсД. В. Многоуровневый подход к разработке алгоритмического и программного обеспечения экзафлопсных суперЭВМ // Вычисл. методы и программирование. 2015. Т. 16. С. 543-556.
11. Мясников Л. В. Мониторинг состояния магматических структур вулкана Эльбрус по наблюдениям литосферных деформаций баксанским лазерным интерферометром: Дис. ... канд. физ.-мат. наук. М., 2012.
12. Гурбанов А. Г., Газеев В. М., Богатиков О. А. и др. Активный вулкан Эльбрус и этапы его геологической истории // Современные методы геолого-геофизического мониторинга природных процессов на территории Кабардино-Балкарии. 2005. С. 94—119.
13. Авдулов М. В., Короновский Н. В. О геологической природе гравитационной аномалии Эльбруса // Изв. АН СССР. Сер. геол. 1962. № 9. С. 67-74.
14. Авдулов М. В. О геологической природе Эльбрусского гравитационного минимума // Вестн. МГУ. Сер. 4. Геология. 1993. № 3. С. 32-39.
15. Собисевич А. Л. Избранные задачи математической геофизики, вулканологии и геоэкологии. Т. 1. М.: ИФЗ РАН, 2012.
Статья поступила 20 августа 2015 г.
Глинский Борис Михайлович, Мартынов Валерий Николаевич, Сапетина Анна Федоровна
Институт вычислительной математики и математической геофизики СО РАН пр. Лаврентьева, 6, Новосибирск, 630090 уптФптБ scc.ru
