CC BY
70
Вычислительная математика
УДК 519.688, 550.34.013.4 DOI: 10.14529/cmse150406
ТЕХНОЛОГИЯ СУПЕРКОМПЬЮТЕРНОГО 3D МОДЕЛИРОВАНИЯ СЕЙСМИЧЕСКИХ ВОЛНОВЫХ
__О __1
ПОЛЕЙ В СЛОЖНО ПОСТРОЕННЫХ СРЕДАХ1
Б.М. Глинский, В.Н. Мартынов, А.Ф. Сапетина
В работе рассматриваются вычислительные технологии решения задач, связанных с моделированием распространения сейсмических волн в неоднородных средах, характерных для вулканических структур, с использованием суперкомпьютерного моделирования в целях создания систем вибросейсмического мониторинга сейсмоопасных объектов. Построена физико-математическая модель магматического вулкана и программная реализация на основе известного численного метода, эффективно использующая архитектуру современного суперкомпьютера, оснащенного GPU. Созданы параллельные 2D и 3D алгоритмы и программы для моделирования распространения упругих волн в сложно построенной среде (2D модель есть сечение исходной 3D модели различными плоскостями и под разными углами) на основе явной конечно-разностной схемы на сдвинутых сетках и метода поглощающих границ CFS-PML. Исследована масштабируемость алгоритмов. Применение разработанной технологии позволяет гораздо эффективней проводить изучение структуры волнового поля, обусловленного геометрией внутренних границ, уточнение его кинематических и динамических характеристик.
Ключевые слова: .мониторинг, 3D моделирование, упругие волны, разностные схемы, гибридный кластер, GPU.
Введение
Необходимость предсказания катастрофических явлений, которые могут быть вызваны готовящейся вспышкой вулканической деятельности, является актуальной задачей. Для решения этой задачи необходимо проведение комплексных и объективных исследований процессов, происходящих на поверхности и внутри вулканической структуры.
Одним из инструментов таких комплексных исследований является активный вибросейсмический мониторинг, при проведении которого, наряду с экспериментальными работами, требуется предварительное и сопутствующее численное моделирование распространения упругих волн в вулканических структурах.
Проведение натурных геофизических экспериментов позволяет получить некоторые представления о скоростных параметрах упругой среды, а также о геометрии изучаемого объекта. В процессе обработки результатов полевого эксперимента часто возникают сложные эффекты, требующие теоретических исследований для их интерпретации и объяснения, которые могут быть выполнены с применением масштабного численного моделирования. Таким образом, создание математических моделей изучаемых объектов и теоретическое моделирование сейсмических полей актуально и помогают в изучении реальных вулканических структур.
1 Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии - 2015».
В связи с тем, что реальный объект исследования имеет сложный рельеф, как правило, не удается поставить площадную систему наблюдения для решения обратной задачи геофизики. Поэтому приходится решать набор прямых задач с целью определения параметров среды, соответствующих экспериментальным наблюдениям на поверхности изучаемых вулканов.
Программный комплекс, решающий поставленную задачу на основе выбранного математического аппарата, должен иметь возможность моделирования упругой неоднородной среды со сложной геометрией. В связи со сложностью и масштабом моделируемой области, решение задачи численного моделирования распространения упругих волн от сосредоточенного источника может требовать значительных вычислительных ресурсов. Поэтому необходима разработка суперкомпьютерных технологий для уменьшения времени расчета и возможностью моделирования «больших» 3D моделей упругих сред. Рассмотрим применение такой технологии на примере моделирования процесса вибросейсмического зондирования вулкана Эльбрус.
Статья организована следующим образом. В разделе 1 проводится построение геофизической модели стратовулкана Эльбрус. Разделы 2 и 3 посвящены постановке и методам решения задачи распространения сейсмических волн в упругой среде. В разделе 4 описан программный комплекс для проведения численного моделирования. Раздел 5 посвящен результатам численного моделирования для тестовой модели, построенной на основе приближенной геофизической модели вулкана Эльбрус из раздела 1. В заключении перечислены основные результаты работы и отмечены планы на будущее.
1. Геофизическая модель стратовулкана Эльбрус
Двуглавый вулкан Эльбрус (абсолютные отметки Западной и Восточной вершин 5642,7 м и 5620 м соответственно) находится на северном склоне Большого Кавказа и располагается на водоразделе рек Малки, Баксана и Кубани, впадающих в Каспийское и Черное моря соответственно. Из-за своего географического положения он как бы «нависает» над плотно заселенными районами Северного Кавказа, прилегающими территориями юга России и севера Грузии. Период образования вулкана 20 млн. лет, последнее извержение 50 год.
Наличие обширного снежно-ледового покрова делает вулкан Эльбрус еще более опасным, так как в случае будущих любой силы и типа извержений, к собственно вулканической опасности, обязательно добавится катастрофическая опасность от образования лахаров и крупномасштабных наводнений, которые могут предваряться валом воды до нескольких десятков метров [1].
На основании данных, приведенных в работах [2-6], а также в работах Н.И. Хитаро-ва (1984), Ю.В. Нечаева (1999) и других авторов построим геофизическую модель стратовулкана Эльбрус со следующими характеристиками. Западная вершина (более старая): высота 5642 м, кратер 500-600 м, глубина 200-300 м; Восточная вершина: высота 5621 м, кратер 200 м, глубина 80 м; между вершинами 3000 м, высота седловины 5322 м; до высоты 4000 м, склоны пологие, но выше 50-600 м, средний наклон 30 0, диаметр у основания 15-18 км, относительная высота 3000 м; вулканическая постройка лежит на гранитном блоке +1 (рис. 1); эффузивные породы, слагают вулканический конус +11; ниже нулевой отметки можно выделить 8 слоев (табл. 1). Зададим верхний магматический очаг в виде эллипсоида с горизонтальными и вертикальными осям 9 и 6 км
(р = 2,1 г/см3; Ур = 2,2 км/с); диаметр бывшего канала 130 м; зададим материнский магматический очаг: эллипсоид с горизонтальными и вертикальными осями 24 и 13 км (р = 1,8 г/см3; Ур = 1,9 км/с), диаметр предполагаемого подпитывающего канала 250 м. Средний канал — цилиндр, диаметром 160 м. Итак, в качестве приближенной модели вулкана Эльбрус можно принять многослойную среду с включениями в виде эллипсов, цилиндров с параметрами, указанными в табл. 1
Рис. 1. Геофизическая модель вулкана Эльбрус и схема вибросейсмического мониторинга
Параметры среды для геофизической модели вулкана Эльбрус
Таблица 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 38 2,85
2. Постановка задачи
Численное моделирование распространения сейсмических волн в сложно построенных упругих неоднородных средах проводится на основе решения полной системы уравнений теории упругости с соответствующими начальными и граничными условиями, записанной в терминах вектора скоростей смещений й = (и, V,W)T и тензора напряжений
] = (]х
^ уу^ т.
,] ,] ,] )т
' ху ^ х^ ут /
В качестве области моделирования рассматривается изотропная 3D неоднородная сложно построенная упругая среда, представляющая собой параллелепипед, одна из граней которого является свободной поверхностью (плоскость z=0).
Рассмотрим прямоугольную декартову систему координат так, что ось Oz направлена вертикально вниз, а оси Ох и Оу лежат на свободной поверхности. Основные уравнения в векторной форме могут быть представлены в следующем виде:
р— = [А]у + Р (,, х, у, т), — = [в ]й
дt
д,
(1)
А =
д_ дх
0 — 0 — 0 —
0 0 д ду д дт 0
д 0 д 0 д
ду дх дт
0 д 0 д д
дт дх ду
в =
(Л + 2ц) —
дх
дх
лд
дх
д
М
М
ду д_ дт 0
лА
ду
д
(л + 2м) —
ду
дА
ду
д
М—
дх
0
дт
лд.
дт
д
(Л + 2м) -
дт
М
дт
М
М
д_
дх д_
ду
где t — время,
р( х, у, т) — плотность,
Л(х, у, ¿) , м(х, у, т) — параметры Ламе.
Предполагается, что параметры упругой среды зависят от трех пространственных переменных X, У и Z.
Начальные условия имеют вид:
]хт 1=0 = 0 , ]ут I, =0 = 0 , ]ху I=0 = 0 , ]хх I=0 = 0 , ]уу 1=0 = 0 , 1=0 = 0 ,
и(х, у, т) 1=0 = 0 , V(х, у, т) I,=0 = 0, П(х, у, т) |,=0 = 0 . (2)
В качестве граничных условий на свободной поверхности выступают:
] I т=0=0 ] | т=0=0 ]„| т=0=0 . (3)
0
0
д
В данной постановке предполагается, что правая часть (массовая сила) имеет следующий вид: Р(,, х, у, т) = Рхг + Руу + Р2к , где г , 7 , к — единичные направляющие векторы координатных осей.
Например, для источника типа «центр давления» получим представление: Р(,,х,у, т) = 5(х - х0)5(у - у0)3(т - т0)/(,), где (х0, у0, т0) координаты источника, а 5( х) — дельта функция.
3. Метод решения задачи
Метод решения поставленной задачи (1)-(3) основан на использовании конечно-разностного метода. Применяется известная конечно-разностная схема, хорошо себя зарекомендовавшая. Алгоритм построения конечно-разностной схемы наиболее подробно изложен в статьях [7-9]. Расчет сеточных коэффициентов в разностной схеме проводится на основе интегральных законов сохранения. Используемая конечно-разностная схема имеет второй порядок аппроксимации по времени и пространству [7], рассматриваются только равномерные сетки. Для примера представим несколько конечно-разностных уравнений используемой схемы:
п+1 п
й , - й
1
п+—
11
п+— п+—
2 2
1/2 2 \ й 1 й 1 пт - п+ - (] 11 -] 1 1 )
Ргу ,к +Рг-\,],к ^г-2(]ххг,},к -]ххЛ,7,к ) . ^7 + ?к ^2к
2 т
1 1
п+— п+—
(] 21 ,1 -] 21
хтг—, 7 ,к +— хтг —, 7 ,к—
2 2 2 2
Ах
■ + -
Ау
+
ат
+ п
^ хг, у ,к '
1 1
п +— п—
(] Л ., 1 -] Л ., 1)
хтг—, 7 ,к-22
п п
й 1 -й 1
w
-М>
хтг—, 7,к— г—, 7,к г—, 7,к-1
= М1 1 , 1 ( 2 . 2— +
г, 7 ,к-
г-1,7,к-
Ах
(4)
(5)
где м1 1 . 1 =
г—, 7,к— 2 2
< 1 ' 4
1 1 1
-+-+-+ -
КМ,7,к 1,7,к ,к-1 1,7,к-1
Для поглощения отражений от границ области моделирования используется вспомогательный метод CFS-PML [10-12]. Его преимущества по сравнению с классическим методом PML заключаются в том, что он дает более «качественную» картину волнового поля для данной задачи, более прост в реализации и экономичен с вычислительной точки зрения.
Поскольку область расчета представляет параллелепипед, а свободная поверхность располагается на верней его грани, то каждая из его границ, за исключением свободной поверхности, окружается поглощающим слоем. Во внутренней области волновое поле рассчитывается по первоначальным конечно-разностным уравнениям, а при попадании волны в зону поглощения происходит расчет по новым формулам с демпфирующими параметрами, описывающими подход к созданию поглощающих границ. Выбор значений демпфирующих параметров для численных расчетов в соответствующих поглощающих слоях проводится на основе результатов работы [11].
2
2
)
г
1
1
4. Программный комплекс для проведения численных расчетов
Для численного моделирования распространения сейсмических волн в моделях сред характерных для вулканических структур предполагается разработать комплекс алгоритмов и программ, позволяющий проводить конструирование модели сложной геометрии и, затем, расчеты одновременно по нескольким программам для этой модели. Он должен включать:
- программу конструирования сеточных моделей сложно построенных сред с включениями, характерными для магматических вулканов;
- программу численного моделирования распространения упругих волн в 3D неоднородных упругих средах с криволинейной свободной поверхностью;
- программу численного моделирования распространения упругих волн в 3D неоднородных упругих средах с прямолинейной свободной поверхностью;
- программу численного моделирования распространения упругих волн в 2D неоднородных упругих средах с криволинейной свободной поверхностью, для заданного сечения рассматриваемой 3D модели;
- программу численного моделирования распространения упругих волн в 2D неоднородных упругих средах с прямолинейной свободной поверхности.
Расчет может производиться одновременно для 3D модели среды и нескольких 2D сечений с целью облегчения интерпретации и получения необходимой информации об особенностях строения волнового поля для выбранной геофизической модели, системы наблюдения и положения источника. В данной реализации рассмотрен вариант алгоритмов и программ для случая прямолинейной свободной поверхности.
Разработка программ велась с учетом особенностей архитектуры гибридного кластера НКС-ЗОТ+GPU ССКЦ ИВМиМГ СО РАН (http://www2.sscc.ru). Он состоит из 40 вычислительных узлов HP SL390s G7, каждый узел содержит два 6-ядерных CPU Xeon X5670 и три карты NVIDIA Tesla M2090 на архитектуре Fermi, у каждой 1 GPU с 512 ядрами и 6 ГБ оперативной памяти GDDR5. Суммарно НКС-ЗОТ+GPU содержит 80 процессоров (480 ядер) CPU и 120 процессоров (61440 ядер) GPU. Пиковая производительность — 85 Тфлопс.
Программы написаны на языке программирования C++ с использованием технологий CUDA (Compute Unified Device Architecture) и MPI (Message Passing Interface).
Разработан построитель моделей среды, который позволяет конструировать на сеточном уровне сложные модели упругих сред на основе идеи Z порядка, описанной в работах [8, 9]. Предполагается, что задана крупноблочная модель среды, составленная из горизонтальных слоев, в вершинах которых задаются параметры среды (Vp, Vs и р). Эти параметры являются непрерывными внутри каждого блока. Разрывы проходят только по граням соседних слоев. Далее происходит интерполяция параметров среды на более «мелкую» расчетную сетку. После того как построена основная сеточная модель трехмерно-неоднородной упругой среды, возможно дальнейшее усложнение ее геометрической структуры. В построенную модель можно «вставлять» различные геометрические объекты со своими упругими параметрами среды, имеющие как численное, так и аналитическое описание.
Для расчетов разработанным программам требуется набор входных файлов с информацией о моделируемой среде и параметрами расчета в определенном формате. Входная информация содержит физические размеры модели, количество шагов дискретизации по трем пространственным направлениям и по времени, информацию о точечном источнике сейсмических волн (его несущая частота и расположение) и информацию о системе наблюдения. Результатом работы программы является набор бинарных файлов, содержащих снимки волнового поля в указанных плоскостях и теоретические сейсмограммы.
4.1. Параллельная реализация и адаптация алгоритма для решения 3D задачи
В ходе работы создана параллельная программа для численного моделирования распространения волн в трехмерных неоднородных упругих средах, реализующая указанный конечно-разностный метод и метод поглощающих слоев и рассчитанная на использование вычислительных узлов гибридного кластера. Подробное описание параллельной реализации и адаптации алгоритмов приведено в статье [13].
Для распараллеливания данной задачи используется декомпозиция области на слои вдоль направления одной из координатных осей (рис. 2). Каждый слой рассчитывается на отдельном узле, где, в свою очередь, он разбивается еще на подслои вдоль другой координатной оси, по числу графических ускорителей на узле.
X
Рис. 2. Декомпозиция расчетной области для гибридного кластера
При такой реализации каждая графическая карта рассчитывает свою сеточную область внутри подслоя на каждом временном шаге независимо от других, за исключением точек, находящихся на границе между двумя соседними областями. Эти точки являются общими для каждой из областей и для продолжения счета необходимо производить обмен информацией об искомых величинах между «соседями». Обмены производятся при помощи технологии MPI. Работа с графическими ускорителями осуществляется с помощью технологии CUDA, в том числе переключение управления между ускорителями. При этом на каждой графической карте параллельная часть кода выполняется как большое количество нитей. Такой гибридный подход обеспечивает высокую степень параллелизма.
Для эффективного использования выбранной для вычислений архитектуры необходима адаптация и оптимизация используемых в ходе моделирования алгоритмов, основанная на знании архитектуры гибридного кластера и его составляющих и необходимых программных средств.
В ходе адаптации выбранных алгоритмов решена проблема выбора сетки блоков и сетки нитей в блоке для осуществления вычислений на GPU. С этим выбором естественно связаны модификации выбранных алгоритмов. Были рассмотрены два наиболее целесообразных (по степени параллелизма) варианта задания сеток и проведен анализ времени исполнения выбранных вариантов. На основе полученных выводов последующая работа велась для полностью трехмерной декомпозиции расчетной области на нити и блоки, размер которых по компоненте x, по возможности должен быть равен или хотя бы кратен длине warp^ (количество физически одновременно исполняющихся нитей на GPU).
Для реализации метода поглощающих границ для каждого поглощающего слоя добавлены свои трехмерные сетки нитей и блоков, соответствующие геометрии каждого поглощающего слоя. Это позволяет избежать увеличения времени счета, значительно сокращая количество ветвлений.
Использование явной конечно-разностной схемы на сдвинутых сетках и разбиение на слои и подслои вдоль разных координатных направлений позволяет обмениваться не всеми компонентами вектора скоростей смещений й = (U,V,W)T и тензора напряжений а = (axx,ayy ,оху ,axz ,ayz )T между каждыми слоем и подслоем, что позволяет сократить количество обменов. Также, если использовать для вычислений только карты, поддерживающие прямое копирование данных с GPU на GPU по PCI-E шине (два GPU из трех на узлах HP SL390s G7), то обмены между вычислительными узлами и графическими картами можно проводить практически параллельно, не используя во втором случае память хоста. Детальная схема обменов представлена на рис. 3.
Рис. 3. Схема обменов данными между узлами и графическими ускорителями
Для уменьшения времени доступа к памяти производится оптимальное расположение используемых трехмерных массивов в глобальной памяти графической карты. Также используется константная память графической карты, в ней хранятся константы, используемые на каждом шаге по времени. Так как решаемая задача оперирует трехмерными массивами, то сложно говорить об использовании разделяемой или текстурной памяти графической карты.
4.2. Преимущества комбинирования решения 2D и 3D задач
Для проведения численного моделирования распространения сейсмических волн в двухмерной упругой среде разработано программное обеспечение на основе программы для трехмерного случая, использующее три графических ускорителя одного вычислительного узла гибридного кластера. В том числе для решения 2D задачи также используется разностная схема на сдвинутых сетках и вспомогательный метод CFS-PML в их двухмерном варианте.
Преимущество данной реализации заключается в возможности провести полномасштабный 2D расчет всего на одном вычислительном узле. Для сравнения, соответствующая 3D задача потребует для расчета решения практически все ресурсы гибридного кластера и намного больше времени для вычислений. Таким образом, предварительные расчеты для характерных 2D сечений исходной 3D модели позволяют за время расчета 3D задачи получить первоначальные результаты. Для примера приведем времена расчета тестовой 3D задачи и соответствующей ей 2D задачи в табл. 2.
Таблица 2
Время расчета тестовых 3D и 2D задач
Задача 3D 2D
Кол-во узлов расчетной сетки (XxYxZxT) 5003x1000 5002x1000
Время расчета 347,2 с 1,49 с
Такой подход значительно ускоряет процесс исследования особенностей волнового поля при моделировании процесса вибросейсмического зондирования вулканической постройки и позволяет провести меньшее количество расчетов для трехмерной задачи, требующих значительно большего числа ресурсов. Комбинирование расчетов 2D и 3D задач позволяет наиболее оптимально использовать архитектуру гибридного кластера.
Построение модели и реальный расчет волнового поля на узле с тремя графическими ускорителями для соответствующих сеток по времени и пространству (6000x9000 узлов по пространству и 25000 шагов по времени) занимает всего 12 минут.
4.3. Исследование времени работы параллельных программ
Для анализа производительности разработанного программного обеспечения исследована его сильная и слабая масштабируемость. Под сильной масштабируемостью будем понимать уменьшение времени счета одного шага одной и той же задачи при использовании большего числа графических ядер в рамках одной графической карты. Ее исследование позволяет понять насколько эффективно применяемые алгоритмы используют архитектуру самой графической карты. Под слабой масштабируемостью будем понимать сохранение времени счета одного шага одного и того же объема задачи при одновременном увеличении количества графических карт.
Слабая масштабируемость исследована для двух вариантов реализации разработанного программного обеспечения для решения 3D задачи, которые направлены на использование распространенных гибридных архитектур, представленных в России. В первом варианте на каждом узле используется две графических карты, поддерживающие прямое копирование данных с GPU на GPU по PCI-E шине, что позволяет, как было сказано выше, практически параллельно производить обмены между вычислитель-
ными узлами и графическими картами. Во втором случае реализации для вычислений используются три карты на каждом вычислительном узле и обмены между GPU проводятся с использованием памяти хоста.
Результаты исследований для программ 3D моделирования приведены на рис. 4, 5 в виде графиков. Из графика на рис. 4 видно, что задача хорошо ложится на архитектуру графической карты (получено ускорение около 40 раз при использовании всех ядер GPU по сравнению с одним ядром GPU).
О 50 100 150 200 250 300 350 400 450 500
Рис. 4. График исследования сильной масштабируемости
Рис. 5. Графики исследования слабой масштабируемости
Из графиков на рис. 5 видно, что эффективность (отношение времени расчета на n узлах в n раз большей задачи к времени расчета на 1 узле исходной задачи) для варианта реализации с тремя GPU на узле падает быстрее, так как процесс обменов не удается распараллелить, в отличии от случая с двумя GPU на узле. Но в целом достигнута эффективность около 92 % при увеличении количества графических карт до 30 штук, что говорит о действенности предложенного процесса обменов между графическими картами и узлами при расчете.
Сравним время, затраченное на численное моделирование 3D модели с использованием для вычислений узлов с GPU, с аналогичным временем для расчета на классическом кластере с CPU. Используем время полномасштабного расчета, проведенного на
вычислительных блэйд-серверах hp ProLiant BL2x220c G5, находящихся в составе НКС-30Т ССКЦ ИВМиМГ СО РАН, приведенное в работе [8]. Параметры расчета указаны в табл. 3. Соответствующее время составило 31 ч 15 мин 17 с (112517 с).
Аналогичное время расчета с теми же размерами расчетной сетки на узлах гибридного кластера с помощью разработанного программного обеспечения составило 2 ч 56 мин (10560 с); полученное ускорение составило 10,66 раза (табл. 4).
Таблица 3
Параметры расчета
по оси X 1677
Кол-во узлов расчетной по оси Y 1059
сетки по оси Z 971
по времени 10313
Таблица 4
Сравнение времени расчетов на гибридном кластере и кластере с классической МРР-архитектурой
Классический MPP кластер Гибридный кластер
Вид узлов кластера HP BL2x220c G5 HP SL390s G7
Количество узлов, использованных для моделирования 20 15
Количество используемых ядер 160 CPU ядер 15360 GPU ядер
Время, затраченное на моделирование 31 ч 15 мин 17 с (112517 с) 2 ч 56 мин (10560 с)
Ускорение в 10,66 раза
5. Результаты численного моделирования
Для иллюстрации работы разработанного комплекса программ численного моделирования распространения сейсмических волн в сложно построенных средах проведены тестовые расчеты, которые продемонстрировали возможности комплекса и проблемы, возникающие при визуализации и интерпретации результатов расчета.
Для расчета в качестве среды взят фрагмент исходной приближенной модели стра-товулкана Эльбрус, описанной в главе 1 и включающий в себя только верхнюю магматическую камеру и прилегающие к ней каналы, расположенные в пятислойной среде. Параметры среды взяты из табл. 1 (слои +1—IV) и описания магматических включений. Система возбуждения состоит из точечного источника типа «центр давления» с частотой 8 Гц, располагающимся вблизи свободной поверхности в левой части расчетной области.
Результаты 3D численного моделирования содержат большой объем информации. В связи с этим, для представления и анализа результатов численного моделирования мы используем теоретические сейсмограммы и мгновенные снимки различных сечений 3D волнового поля плоскостями, проходящие через линии возможной системы наблюдений.
На рис. 6 представлены снимки волнового поля в плоскости, проходящей через точечный источник и ось симметрии верхнего канала.
А1 (Т = 1,61 с) В1 (Т = 1,61 с)
А2 (Т = 3,22 с) В2 (Т = 3,22 с)
А3 (Т = 4,83 с)
В3 (Т = 4,83 с)
Рис. 6. Мгновенные снимки сечения 3D волнового поля в различные моменты времени, компонента и, плоскость XZ. А1-А4 снимки полученные в ходе 3D моделирования, В1-В4 снимки полученные в ходе 2D моделирования
Слева буквами А1-А4 обозначены мгновенные снимки, полученные в ходе решения 3D задачи, а справа буквами В1-В4 обозначены снимки, полученные в ходе 2D моделирования. Визуализация снимков проведена с помощью программы Aspis, разработанной
03721971
в ОАО «Сибнефтегеофизика». Фактически, данные снимки являются иллюстрацией нашего предположения о том, что 2D моделирование может быть использовано для возможного построения 3D моделей среды.
Из приведенных мгновенных снимков видно, что волновое поле имеет сложную картину и существенно зависит от геометрии, размеров и свойств эллиптических включений. Для проведения геофизической интерпретации полученных теоретических сейсмограмм, даже для этого фрагмента исходной приближенной модели вулкана Эльбрус, потребуется проведение большой серии вычислительных экспериментов. Для слоистой среды, содержащей включения, у которых может меняться ориентация и форма их границ (например, присутствие в модели пересечения нескольких объектов различной природы под различными углами), обусловленные приближением к реальной модели среды, задача значительно усложняется. Тем не менее, представленные эксперименты показывают, что численное моделирование может дать значительную информацию для организации проведения натурных экспериментов и интерпретации результатов экспериментальных наблюдений, полученных в ходе вибросейсмического мониторинга.
Заключение
В работе предложена суперкомпьютерная технология решения задачи распространения упругих волн в сложно построенных средах, учитывающая различные по размерности численные методы исследования неоднородной среды и архитектуру суперкомпьютера.
Предложен комплекс параллельных алгоритмов и программ для гибридного кластера с графическими ускорителями, реализующий распространения упругих волн в 2D и 3D моделях в сложно построенных средах. Причем, фактически созданы две программы, позволяющие проводить расчеты на гибридных кластерах как с двумя, так и стремя графическими картами на узле.
Проведенные численные эксперименты для упрощенной модели вулкана Эльбрус показывают, что волновое поле от точечного источника имеет сложную картину и существенно зависит от геометрии, размеров и свойств эллиптических включений. Особенности распространения упругих волн в средах такого типа могут быть использованы при вибросейсмическом мониторинге вулканических структур. Проведение серии расчетов позволяет выбрать модель вулканической структуры адекватную результатам натурных наблюдений по кинематическим и динамическим характеристикам.
В дальнейшем предполагаются развитие предложенной суперкомпьютерной технологии численных расчетов и последующее проведение численных экспериментов для моделей с криволинейной свободной поверхностью.
Работа поддержана проектом РФФИ №13-07-00589. Литература
1. Лаверов, Н.П. Новейший и современный вулканизм на территории России. / Н.П. Лаверов [и др.] — М.: Наука, 2005. — 604 с.
2. Мясников, A.В. Мониторинг состояния магматических структур вулкана Эльбрус по наблюдениям литосферных деформаций баксанским лазерным интерферометром [Текст] : дис. / А.В. Мясников — Москва, 2012.
3. Гурбанов, А.Г. Активный вулкан Эльбрус и этапы его геологической истории / А.Г. Гурбанов, В.М. Газеев, О.А. Богатиков [и др.] // Современные методы геолого-геофизического мониторинга природных процессов на территории Кабардино-Балкарии. — 2005. — С. 94-119.
4. Авдулов, М.В. О геологической природе гравитационной аномалии Эльбруса / М.В. Авдулов, Н.В Короновский // Изв. АН СССР. Сер.геолог. — 1962. — № 9. — С. 67-74.
5. Авдулов, М.В. О геологической природе Эльбрусского гравитационного минимума / М.В. Авдулов // Вестник МГУ. Сер. 4. Геология. — 1993. — № 3. — С. 32-39.
6. Собисевич, А.Л. Избранные задачи математической геофизики, вулканологии и геоэкологии. / А.Л. Собисевич — М.: ИФЗ РАН, 2012. — T. 1. — 512 с.
7. Bihn, M. Stable Discretization Scheme for the Simulation of Elastic Waves / M. Bihn, T.A. Weiland // Proceedings of the 15th IMACS World Congress on Scientific Computation, Modelling and Applied Mathematics (IMACS 1997). — 1997. — Vol. 2. — P. 7580.
8. Глинский, Б.М. Численное моделирование и экспериментальные исследования грязевого вулкана «Гора Карабетова» вибросейсмическими методами / Б.М. Глинский, Д.А. Караваев, В.В. Ковалевский, В.Н. Мартынов // Вычислительные методы и программирование. — 2010. — Т. 11. — С. 95-104.
9. Караваев, Д.А. Параллельная реализация метода численного моделирования волновых полей в трехмерных моделях неоднородных сред / Д.А. Караваев // Вестник Нижегородского университета им. Н. И. Лобачевского. — 2009. — № 6(1). — С. 203-209.
10. Drossaert, F. Complex frequency shifted convolution PML for FDTD modeling of elastic waves / F. Drossaert, A. Giannopoulos // Wave Motion. — 2007. — Vol. 44. — P. 593604. DOI: 10.1016/j.wavemoti.2007.03.003.
11. Komatitsch, D. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation / D. Komatitsch, R. Martin // GEOPHYSICS. — 2008. — Vol. 73, No. 4. — P. T51-T61. DOI: 10.1190/1.2939484.
12. Hastings, F.D. Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation / F.D. Hastings, J.B. Schneider, S.L. Broschat // J. Acoust. Soc. Am. — 1996. — Vol. 100(5). — P. 3061-3069. DOI: 10.1121/1.417118.
13. Сапетина, А.Ф. Численное моделирование распространения сейсмических волн в сложно построенных средах на гибридном кластере / А.Ф. Сапетина // Проблемы прочности и пластичности. — 2014. — T. 76, № 4. — С. 288-296.
Глинский Борис Михайлович, д.т.н., заведующий лабораторией, Сибирский суперкомпьютерный центр Института вычислительной математики и математической геофизики СО РАН, заведующий кафедрой вычислительных систем, Новосибирский государственный университет (Новосибирск, Российская Федерация), gbm@opg.sscc.ru.
Мартынов Валерий Николаевич, с.н.с. лаборатории численного моделирования сейсмических полей, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), vnm@nmsf.sscc.ru.
Сапетина Анна Федоровна, аспирант, инженер лаборатории сибирского суперкомпьютерного центра, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), afsapetina@gmail.com.
Поступила в редакцию 15 апреля 2015 г.
Bulletin of the South Ural State University Series "Computational Mathematics and Software Engineering"
2015, vol. 4, no. 4, pp. 101-116
DOI: 10.14529/cmse150406
TECHNOLOGY OF SUPERCOMPUTER SIMULATION OF SEISMIC WAVE FIELDS IN COMPLICATED MEDIA
B.M. Glinskiy, Institute of Computational Mathematics and Mathematical Geophysics of Siberian Branch of Russian Academy of Sciences (Novosibirsk, Russian Federation) gbm@opg.sscc.ru,
V.N. Martynov, Institute of Computational Mathematics and Mathematical Geophysics of Siberian Branch of Russian Academy of Sciences (Novosibirsk, Russian Federation) vnm@nmsf. sscc.ru,
A.F. Sapetina, Institute of Computational Mathematics and Mathematical Geophysics of Siberian Branch of Russian Academy of Sciences (Novosibirsk, Russian Federation) afsapetina@gmail.com
The paper considered computing technology solving problems related to the modeling of seismic wave propagation in inhomogeneous media typical of volcanic structures using supercomputer simulations in order to create systems of vibroseis monitoring for quake-prone objects. The physico-mathematical model of the magmatic volcano is constructed and software implementation on the basis of the known numerical method that effectively using the architecture of modern supercomputers equipped with GPU is developed. The parallel 2D and 3D algorithms 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) on basis of the explicit finite-difference scheme for the shifted grids and CFS-PML method of absorbing boundaries is developed. Scalability of algorithms is investigated. The application of the developed technology allows for much more efficient to carry out studies of the structure of the wave field due to the geometry of the internal boundaries and refinement of its kinematic and dynamic characteristics.
Keywords: monitoring, 3D simulation, elastic waves, finite difference schemes, hybrid cluster, GPU.
References
1. Laverov N.P., et al. Noveyshiy i sovremennyy vulkanizm na territorii Rossii [Modern and Holocene volcanism in Russia]. Moscow, Science, 2005. 604 p.
2. Myasnikov A.V. Monitoring sostoyaniya magmaticheskikh struktur vulkana El'brus po nablyudeniyam litosfernykh deformatsiy baksanskim lazernym interferometrom [Monitoring the State of the Magmatic Structures of Elbrus Volcano Based on Observation of
Lithosphere Strains by observation of the Baksan Laser Interfrometer]. dissert. Moscow, 2012.
3. Gurbanov A.G., Gazeev V.M., Bogatikov O.A. and etc. Aktivnyy vulkan El'brus i etapy ego geologicheskoy istorii [Active volcano Elbrus and the stages of its geological history] // Modern methods of geological and geophysical monitoring of natural processes on the territory of Kabardino-Balkaria. 2005. P. 94-119.
4. Avdulov M.V., Koronovskiy N.V. O geologicheskoy prirode gravitatsionnoy anomalii El'brus [About the biological nature of Elbrus gravity anomaly] // The news of AS USSR. Geology Series. 1962. No. 9. P. 67-74.
5. Avdulov M.V. O geologicheskoy prirode El'brusskogo gravitatsionnogo minimuma [About the biological nature of the Elbrus gravity minimum] // Moscow University Geology Bulletin. 1993. No. 3. P. 32-39.
6. Sobisevich A.L. Izbrannye zadachi matematicheskoy geofiziki, vulkanologii i geoekologii [Selected problems of Mathematical Geophysics and Volcanology and Geoecology]. Moscow, IPE RAS, 2012. Vol. 1. 512 p.
7. Bihn M., Weiland T.A. Stable Discretization Scheme for the Simulation of Elastic Waves // Proceedings of the 15th IMACS World Congress on Scientific Computation, Modelling and Applied Mathematics (IMACS 1997). 1997. Vol. 2. P. 75-80.
8. Glinskiy B.M., Karavaev D.A., Kovalevskiy V.V., Martynov V.N. Chislennoe modeliro-vanie i eksperimental'nye issledovaniya gryazevogo vulkana «Gora Karabetova» vi-broseysmicheskimi metodami [Numerical modeling and experimental research of the «Karabetov Mountain» mud volcano by vibroseismic methods] // Numerical Methods and Programming. 2010. Vol. 11. P. 95-104.
9. Karavaev D.A. Parallel'naya realizatsiya metoda chislennogo modelirovaniya volnovykh poley v trekhmernykh modelyakh neodnorodnykh sred [Parallel implementation of wave field numerical modeling method in 3D models of inhomogeneous media] // Vestnik of Lobachevsky State University of Nizhni Novgorod. 2009. No. 6(1). P. 203-209.
10. Drossaert F., Giannopoulos A. Complex frequency shifted convolution PML for FDTD modeling of elastic waves // Wave Motion. 2007. Vol. 44. P. 593-604. DOI: 10.1016/j.wavemoti.2007.03.003.
11. Komatitsch D., Martin R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation // GEOPHYSICS. 2008. Vol. 73, No. 4. P. T51-T61. DOI: 10.1190/1.2939484.
12. 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. Vol. 100(5). P. 3061-3069. DOI: 10.1121/1.417118.
13. Sapetina A.F. Chislennoe modelirovanie rasprostraneniya seysmicheskikh voln v slozhno postroennykh sredakh na gibridnom klastere [Numerical simulation of seismic wave propagation in a complicated medium on the hybrid cluster] // Problems of Strength and Plasticity. 2014. Vol. 76, No. 4. P. 288-296.
Received April 15, 2015.
