Моделирование упругих волн в средах со сложной топографией свободной поверхности Текст научной статьи по специальности «Математика»

УДК 519.6

DOI 10.25205/1818-7900-2018-16-4-153-166

П. А. Титов

Институт вычислительной математики и математической геофизики СО РАН пр. Академика Лаврентьева, 6, Новосибирск, 630090, Россия

[email protected]

МОДЕЛИРОВАНИЕ УПРУГИХ ВОЛН В СРЕДАХ СО СЛОЖНОЙ ТОПОГРАФИЕЙ СВОБОДНОЙ ПОВЕРХНОСТИ *

Моделирование - неотъемлемая часть исследования процессов распространения сейсмических волн в различных средах. Широко используемый способ разбиения области на прямоугольные ячейки обладает недостатком: при сложном строении свободной поверхности (например, гора) возникают эффекты при отражении волны от этой поверхности, связанные с тем, что граница аппроксимируется ступенчатой функцией для численного решения задачи. В данной работе предложен иной подход к дискретизации области: построение криволинейной сетки, хорошо согласующейся с геометрией свободной поверхности. Предложен алгоритм численного моделирования на основе численного решения линейной 2Б-системы теории упругости, записанной в смещениях, с использованием криволинейной сетки и пошагового метода Лагерра. Представлены результаты моделирования. Также были реализованы две параллельные версии алгоритма, проведены расчеты на различных многоядерных системах (на классической многопроцессорной архитектуре, а также архитектуре с использованием сопроцессоров Intel Xeon Phi - "РСК ПетаСтрим"). Представлены сравнительные тесты ускорений на разных архитектурах, а также описаны особенности распараллеливания алгоритма.

Ключевые слова: теория упругости, упругие волны, криволинейная поверхность, криволинейные сетки, метод Лагерра, Intel Xeon Phi, PetaStream.

Введение

Эффективным инструментом исследования процессов распространения волн в различных моделях сложно построенных сред является математическое моделирование. Одним из известных методов численного моделирования распространения упругих волн в случае сложно построенных 2Б-сред является разностный метод. При использовании данного метода неизменно сталкиваются с проблемой возникновения дифракционных волн при отражении волны от свободной поверхности. Связано это с тем, что для применения данного метода необходимо разбиение физической области на квадратные или прямоугольные ячейки, вследствие чего свободная поверхность заменяется ступенчатой функцией. Пример использования разностного метода приведен на рис. 1 и 2.

Подобные дифракционные волны вносят погрешность в решение. Таким образом, была поставлена задача об избавлении от дифракционных волн. Для этого в работе предлагается

* Работа выполнена в рамках государственного задания ИВМиМГ СО РАН (проект № 0315-2016-0009), а также при поддержке грантов РФФИ № 16-07-00434 и № 16-01-00455 А.

Титов П. А. Моделирование упругих волн в средах со сложной топографией свободной поверхности // Вестн. НГУ. Серия: Информационные технологии. 2018. Т. 16, № 4. С. 153-166.

ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2018. Том 16, № 4 © П. А. Титов, 2018

Рис. 2. При отражении от свободной поверхности видны дифракционные волны, идущие вслед за отраженными волнами (возле поверхности) (Рисунок предоставлен Д. А. Караваевым, ИВМиМГ СО РАН)

использовать криволинейные сетки. В последнее время все большую популярность приобретает моделирование с использованием криволинейных сеток. Это гибкий инструментарии, нашедший применение в численном решении сингулярно-возмущенных задач, моделировании потоков плазмы в камере Токомака и др. Способы построения сеток, а также их применение для решения упомянутых задач приведены в [1].

Данная работа особенно важна для дальнейшей разработки 3Б-алгоритма, нацеленного на развитие методов исследования областей, характерной особенностью которых является сложное строение топографии свободной поверхности. К таковым можно отнести магматические вулканы, а также территории арктического шельфа и Сибири (Баженовская свита), богатые залежами углеводородов. Существующие методы работы с такого рода задачами не позволяют получить желаемую точность моделирования ввиду большой трудозатратности и особенностей применяемых методов, таких как метод конечных элементов (МКЭ).

Построение криволинейной сетки

Существует много методов построения криволинейной сетки. Нкоторые из них описаны в работах [1-12]. Криволинейная сетка в «физической» области получается в результате взаимно-однозначного отображения равномерной сетки «расчетной» области. Таким образом, необходимо установить взаимно-однозначное соответствие между «физической» и «расчетной» областью. «Физическая» область находится в пространстве переменных (х,у),

а «расчетная» - в пространстве переменных (д, г). Для этого был использован метод

трансфинитной интерполяции [1, с. 53]. Автором были подобраны управляющие функции, позволяющие определять свойства криволинейной сетки.

В данной работе важным моментом построения сетки является ее ортогональность к свободной поверхности, что накладывает ограничение на форму свободных поверхностей, с которыми можно работать при помощи данного метода построения сетки, а именно: поверхность должна быть гладкой функцией у = G (х) первого порядка на промежутке, на котором она построена для данной области, а также для любых двух разных у1 и у2 не существует х, такого что у1 = G (х) и у2 = G (х). Пример на рис. 3.

Рис. 3. Пример областей

Таким образом, задачу, поставленную в переменных (х,у), можно переписать в переменных (д, г ).

Постановка и преобразование задачи

Моделирование проводится на основе численного решения полной линейной системы теории упругости, записанной в смещениях. В «физической» области в переменных (х, у) линейная система уравнений в смещениях принимает вид

д2и д (, . дм Л д

(2Ц + А,)— + Х— + — дх ду) ду

2

дх

( ( Ц

V V

дм ди — + —

дх ду

д2 м

д( дх

( ( Ц

V V

дм ди — + —

дх ду

)

ЛЛ

)

д (\ .лди дм Л — (2ц + А,)— + Х—

ду V дх ду

+ Е. = 0,

+ Еу =0.

Р

Граничные условия на дГ и начальные условия:

I I п дw ды

ы дг= ы дг= 0, = =

дt д/

(2)

Условия на свободной поверхности д5":

Л

|(2ц + Л) +

^ ' дх ду

( (ды ды ^ ( Ц

( ( ды ды ^

Ц

V V

дх ду

= 0,

//

дх ду

//

. ч ды ды (2ц + Я)—+ Я—

дх ду

\

(3)

= 0,

где ы(х,у), ы(х,у) - смещения вдоль координатных осей ОХ и ОУ; Х(х,у), ц(х,у) -параметры Ламэ; р(х,у) - плотность; дГ - граница области, не включающей свободную поверхность; дS - свободная поверхность; [пх,пу] - единичный вектор нормали к свободой поверхности дS; , р - правая часть, отвечает за работу источника возмущений

в среде. В данной работе используется источник типа «центр давления», работает в течение 2-х секунд.

К =

Ру =

д5(х дху у°) (п( л (/ -1)) + 0,88т(2л (/ -1)) + 0,28т(3л(/ -1))), 0

0, 1 > 2,

д5(х ду у у0) (п ( л (/ - 1)) + 0,8 8т (2 л (/ -1)) + 0,2 8Ш (3 л (/ -1))), 0

< 1 < 2,

< 1 < 2,

0, 1 > 2.

Аналогично работе [13] после преобразования системы (1)-(3) получаем задачу вида

•¿Р 2

д/2 дд

ды д [(2Ц+^)(ч*дд + гхдг)ы+*(дч + гудг)ы] + 3ду ц(( дд + Гх дг ) +(ду дд + Гудг )ы )

+-

дг +3гу

¿Гх [(2Ц + Я)(дд +Гхдг )ы + ^Чудд +Гудг )]-

г д2 ы д

3 Р—т = —

д/2 дд

Ц((дд + гхдг ) + (дудд + гудг )ы)]] + Рх = В (ы, + К: 3Чх [Ц((дд + гхдг ) ы + (дудд + гудг )ы)] +

+3Чу [(2Ц + ^)( Чх дд +гх дг) ы + х( ду дд +гу дг) ы]

(4)

+

дг

3г

ц((ч д + гд ) + (ч д + гд )ы ) 1 +

~\д д х г' V у д у г' /]

+ 3гу [(2Ц+^)( Чх д д + гхд г ) ы + ^ ( Чу д д + гу д г ) ™ ]] + Ру = В2 (ы , ) +

Начальные и граничные условия:

ы 0Г = ы 0Г = 0,

ды ды д/ д/

= 0.

(5)

д

д

На свободной границе:

= 0,

гх [(2Ц + Ь)(хдд + гхд)и + х(дудд + Гудг )м] +

+Гу [ц((Ухдд + Гхдг )м + (Чудд + Гудг )и)

ц((хдч + гхдг ) + (Чудч + гудг )и) +

+Г [(2Ц + Ь)( Чх д Ч + гх д г) и + х( Чу д Ч + гуд г) = 0

(6)

где г =

г

гу =

(г

Чх = —, гх =

У х

А

J

Чу =■

J

Гу J

J = ч г — ч г -

1х у 1у х

Якобиан проебразования.

Пошаговый метод преобразований Лагерра по времени

Впервые данный метод был применен в работе [14]. В отличие от классических преобразований Фурье и Лапласа применение преобразования Лагерра приводит к решению задачи для системы дифференциальных уравнений, не зависящей от параметра разделения, роль которого в данном случае играет номер проекции Лагерра. В результате применения этого преобразования получается задача для проекций Лагерра с правой частью, ре-куррентно зависящей от проекций меньших номеров. После применения конечно-разностных аппроксимаций по пространственным переменным, решение исходной задачи сводится к решению системы линейных алгебраических уравнений со многими правыми частями, для решения которых можно использовать различные современные подходы к решению линейных систем.

В работе рассматривается модификация метода, а именно преобразование Лагерра используется на последовательности конечных интервалов по времени. Полученное решение в конце одного временного отрезка используется в качестве начальных данных для решения задачи на следующем временном отрезке и т. д. Связано это с тем, что для получения решения с хорошей точностью через длительный промежуток времени необходимо очень большое количество проекций Лагерра, что критично скажется во время численной реализации. Подробно пошаговый метод преобразования Лагерра описан в работе [15].

При использовании данного подхода появляется необходимость выбора 4-х параметров: количество проекций Лагерра N масштабного множителя, необходимого для аппроксимации решения функциями Лагерра к, экспоненциального коэффициента весовой функции ц, использующейся для нахождения решения на конечном временном интервале, и длительности этого интервала т. Способ выбора этих параметров предложен в работе [15]. Мы же ограничимся лишь основными терминами и непосредственно видом задачи после примененного к ней преобразования Лагерра.

Полиномы Лагерра имеют вид

Ь„

¿о (() = 1,

Ц( () = 1 — (, () = 1+1 ( +1 —') (( — кЬ*—1 ((), к > 1, 4 (()=—1£1, (), к >1.

(7)

(8)

(9) (10)

г

х

г

х

У

х

г

Функции Лагерра

¡п ^) = в~2Ьп ^), п > 0,

(11)

образуют полную ортонормированную систему в Ь0 (0,да). Преобразование

да

^п =1V ( к

(12)

называется прямым преобразованием Лагерра, а уп - проекциями Лагерра функции V (t) . Обратное преобразование Лагерра:

да

V() = ЕVn¡n (ht). (13)

п=0

Идея метода состоит в получении приближенного значения

N

V(т) - X Vn¡n (^) (14)

п=0

и значения

V (т) = (т), V (т) = ци (т) + Л> (т)

(15)

выбираем в качестве начальных данных для решения задачи при t >т. Иначе говоря, теперь вместо нахождения неизвестной функции V (t) мы ищем значения проекций Лагерра для этой функции.

Применение преобразования Лагерра к задаче (3)-(6)

Учитывая уравнения (7)-(15), задача (3)-(6) сводится к решению линейной системы уравнений:

V+—| +А 2

V+—| +А 2

V+—| +А 2

(Ь2 . ^

— + 2vh 2

V0 + hv1 + /о:

v1 = -(2 + 2уЬ )0 + ^ 3 Ь2 + 2vh ^ V0 + ЬУ1 + /1,

(п - 2^-1 + Уп-2 ) + 2УЬ (-1 - Уп-2 ) + ЬЧ-1 = / - 2/п-1 + /,

п > 2.

(16)

На дГ:

Уп| <г= 0, п = 0,..., N.

На дБ:

К [(2^)( д„ +Г дг ) +^(Чу дЧ +Гу дг ) ] +

(17)

(18)

= 0,

Ч [Юд*д* + Гхдг ) +(дУд* + ГУдг )Ып )

ц((х д * + Г д г ) + (*у д * + Гу д г ) )] +

+гу [(2ц + Яхд* + гхдг)ип + х(дуд* + гудг) ып ] = 0.

В (16) вместо Уп подставляются ип, ып - п-е проекции Лагерра для и и ы соответственно.

В1 В2

Вместо оператора л подставляются — для ип и — для ып соответственно. N - число

Jp Jp

проекций Лагерра, V, к- параметры преобразования Лагерра. Таким образом вместо изначальной системы из 2-х уравнений, записанной в смещениях, получаем систему из 2 ( +1) уравнений, записанную в проекциях Лагерра.

Далее для полученной системы строится разностный аналог на равномерной сетке размером М х К ячеек. Разностная схема для В, (и,ы), В2 (и,ы) описана в работе [13,

уравнения (3.3), (3.4)], за исключением дискретных граничных условий на свободной поверхности дБ и на дГ, которые предложены автором данной работы. Упростив (18), получим на дБ:

а

¡к

1 дд

а

дип

дг

а

дып д*

ды

а

„ ди „ ди

Р^ + Р2"Тп + Р3 д

д* дг д*

дып

дг ды

_п

дг

= 0,

= 0.

(19)

Из (19) получаем разностный аналог:

„ = и*1 - V2

р1а4 ~а1р4 Ча2Р4 -Р2а4 У

(2 (ип+1Д - ип-1,1 )- V2 (( - ип-,2 )

,1/2

-1/2

^Р3а4 -а3Р4 ^

ы.

= ыц -1/2

(,0 п 1 '

^р1а2 -а1р2 ^ Ча4Р2 -Р4а2 У

,1/2

Ча2р4 -р2а4 У ((2 ((+1Д - ип-1Д )- V2 ( +1,2 - ип,-1,2 ))

(2 (ып+1Д - ып -1Д )-V2 (ып +1,2-ып -12 ))

,1/2

-1/2

^р3а2 а3Р2 ^ Ча4р2 -р4а2 У

(2(ып+1Д-1Д ))(ып +1,2)).

Д/2

Здесь используются так называемые «фиктивные» узлы. На дГ из (17) получаем разностный аналог:

ип = ип = ып = ы„ = 0, ип = ы„ = 0, п = 0,..., N.

п.,] пм,, 1,] пм,, ' п,к п,к ' '

В итоге получаем разностный аналог системы (16)-(18) со вторым порядком аппроксимации по пространственным и временной переменным.

Для численного решения был задействован метод простой итерации. Подробно суть метода изложена в [16].

Выбор этого метода обусловлен тем, что система для численного решения из 2 (+1)

уравнений, переписанная в матричном виде, будет обладать матрицей с большим диагональным преобладанием, что обеспечит быструю сходимость метода.

Далее, на примере V», (вместо него подставляются и wn ) опишем общий алгоритм

действии.

0. На интервале [0, т] начальные данные выглядят так: У. 0 = 0), у1.. = У . (0).

0) и У Д 0) берутся из начальных данных задачи.

1. Получив все значения V» при помощи метода простой итерации, находим

(т) = Т, N=0 {тН ).

2. Отсюда, используя формулы V= с= ^^=0сп1п (Н) с0 = Н- V(0^, сп =

Н

= сп-1 + ^ (» + vn-l), п = - •, найдем V..

3. Затем у.(т) = в^и(т), У_1. (т) = ^. (т) + e^тVi. (т) и, наконец, значения У0 = Уи(т), У]. = У^. (т) берутся в качестве начальных данных для следующего временного интервала.

Присваиваем V0. = у.°, V1,. =-у +у1..

4. vi j (х) и (х) будут начальными данными для следующего временного интервала. Далее, аналогично находим v^ (2х), V^ (2х) и т. д.

Программная параллельная реализация алгоритма

При помощи языка Fortran и пакета MPI было создано 2 варианта параллельной программы и проведены численные расчеты на кластере ССКЦ НКС-30Т (на узле G7 2 процессора Intel Xeon E5670 и 24 ГБ оперативной памяти. http://www.sscc.icmmg. nsc.ru/hardware.html) с классической многопроцессорной архитектурой, а также на кластере МВС-10П МП с архитектурой РСК ПетаСтрим (8 модулей по 8 сопроцессоров Intel Xeon Phi 7120D, 61 ядеро, 244 потока, RAM 16 Gb DDR5. На каждом модуле также стоит процессор Intel Xeon E5-2667 v2. Твердотельные накопители Intel SSD DC S3500. http:// www.jscc.ru/resources/hpc/). Декомпозиция области для реализации на CPU и на Intel Xeon Phi представлена на рис. 4.

Особенностью архитектуры РСК ПетаСтрим является то, что пользователю доступны для работы только сопроцессоры Intel Xeon Phi. Между собой они обмениваются при помощи MPI-пересылок, а внутри каждого Phi также осуществляется распараллеливание по потокам при помощи OpenMP. Подробнее с особенностями программирования для Intel Xeon Phi можно ознакомиться в работах [17-21]. Далее будут приведены сравнительные тесты ускорений на примере конкретной модели среды.

Для решения системы линейных уравнений используется метод простых итераций, что позволяет «прореживать» обмены между процессами, так как это не влияет на сам факт сходимости итерационного процесса, а только на ее скорость. В ходе расчетов для модели, описаной далее, выяснилось, что относительно временных затрат на расчеты выгоднее сделать суммарно вдвое большее число итераций и вдвое меньшее число обменов. Другими словами, соседние процессы обмениваются раз в четыре итерации. Это ускорило работу программы в среднем в 1,5 раза, при этом точность решения пострадала незначительно -в тестовых расчетах расхождение решений по амплитуде, полученных с прореживанием обменов и без прореживания, составило порядка 0,001 %.

Математическая модель

Область - однородная среда с горой на поверхности (рис. 5).

Рис. 4. Декомпозиция расчетной области для реализации на CPU (а) и Intel Xeon Phi (б)

Параметры среды : р = 1, X = 0,5, ц = 0,25. На рис. 6 размеры области даны в километрах. Соответствующая ей «расчетная» область - прямоугольник размером 64 х 32 км, дискретизация области - 5 120 х 2 060 ячеек. Источник типа «центр давления» работает в течение 2-х секунд. В уравнениях он представлен двумя компонентами (р, р ):

F =

F =

Ô5(ô ^y Уо) (sin( л (t -1)) + 0,8sin(2л (t -1)) + 0,2sin(3 л(t -1))), 0 < t

0, t > 2,

< t < 2,

Ô5(x Ôy У У0 ) ( (л (t -1)) + 0,8sin (2л (t -1)) + 0,2 sin (3 л (t -1))), 0 < t

< t < 2,

0, t > 2.

а

б

Рис. 5. Область с «горкой». Точка - расположение источника

Рис. 7. Снимки волнового поля в последовательные моменты времени 2, 7, 12, 17, 22, 27, 32, 37 с

Важным свойством криволинейной сетки яляется ее ортогональность к свободой поверхности (рис. 6). Следует также упомянуть, что для двух разных форм поверхности не всегда можно использовать один и тот же алгоритм построения сетки. Зачастую возникает необходимость подбирать параметры управляющих функций занова, чтобы не допустить «перехлеста» ячеек криволинейной сетки. Аналитический метод построения позволяет делать это быстро. Программный генератор сетки написан на языке Fortran и также предложен автором.

Результаты численного моделирования

На рис. 7 можно наблюдать отсутствие дифракций при отражении волны от поверхности, что подтверждает эффективность подхода, предложенного в данной работе.

Результаты тестов

Далее приведены сравнительные тесты для разных архитектур (рис. 8). Модель описана в разделе «Математическая модель». Тесты для CPU выполнены на сервере G7 (2 х 6-ядерных Intel Xeon E5670 и 24 Гб оперативной памяти на узле) кластера ССКЦ НКС-30Т, при этом технология hyperthreading отключена. Для Intel Xeon Phi - на кластере МВС-10П. В случае Xeon Phi не проводилось никаких дополнительных процедур оптимизации вроде выравнивания массивов или изменения числа потоков на ядро (по умолчанию их 4). В тесте на масштабируемость алгоритма на каждый сопроцессор Phi приходился объем вычислений, соответствующий нашей модели (5 120 х 2 060 ячеек).

Рис. 8. Тесты. Слева - для CPU, справа - для PetaStream

Как видно из результатов, даже один узел с 12 ядрами CPU (на 1 ядро - 1 MPI-поток) справляется быстрее, чем один Intel Xeon Phi. При этом заметно, что в случае PetaStream алгоритм быстро «уходит в насыщение», что может свидетельствовать о недостаточной оптимизации MPI-обменов между соседними Phi на уровне программного обеспечения кластера.

Заключение

В работе представлен алгоритм и результаты численного моделирования с использованием этого алгоритма. Также непосредственно автором создано 2 параллельных реализации алгоритма, программно реализован генератор криволинейной сетки. Проведены сравнительные тесты для разных архитектур на примере конкретной модели.

Новизна работы заключаются в следующем: предложен новый способ аппроксимации граничных условий, впервые моделирование сейсмических полей проводилось с использованием локально-ортогональной криволинейной сетки, построенной аналитическим способом.

Проведенные тесты показали отсутствие дифракции при отражении от криволинейной свободной поверхности. В дальнейшем эта работа послужит основой для полноценной 3D-задачи, результаты которой могут быть применены для решения реальных задач.

Автор выражает благодарность В. Н. Мартынову за консультации в ходе выполнения работы, а также ССКЦ РАН и МСЦ РАН за предоставленное оборудование для численных расчетов.

Список литературы

1. Лисейкин В. Д., Рынков А. Д., Кофанов А. В. Технология адаптивных сеток для решения прикладных задач: Монография / Новосиб. гос. ун-т. Новосибирск, 2011. 160 с. ISBN 978-594356-981-4

2. Хакимзянов Г. С., Шокин Ю. И. Лекции по разностным схемам на подвижных сетках. Часть 2. Задачи для уравнений в частных производных с двумя пространственными переменными. Казань, 2006.

3. Васева И. А., Кофанов А. В., Лисейкин В. Д., Лиханова Ю. В., Харитончик А. М. Технология построения пространственных адаптивных сеток // Журнал вычислительной математики и математической физики. 2010. Т. 50, № 1. С. 1-19.

4. Васева И. А., Лисейкин В. Д. Применение метода конечных элементов для построения адаптивных сеток // Вычислительные технологии. 2011. Т. 16, № 6. С. 30-35

5. Глассер А., Лисейкин В. Д., Китаева И. А. Контролирование свойств разностных сеток с помощью мониторной метрики // Журнал вычислительной математики и математической физики. 2005. Т. 45, № 8. С. 1466-1483.

6. Глассер А., Лисейкин В. Д., Шокин Ю. И., Васева И. А., Лиханова Ю. В. Построение разностных сеток с помощью уравнений Бельтрами и диффузии. Новосибирск, Наука, 2006.

7. Лисейкин В. Д., Мороков Ю. И., Васева И. А., Лиханова Ю. В. Некоторые аспекты построения адаптивных сеток // Журнал вычислительной математики и математической физики. 2008. Т. 48, № 9. С. 1638-1658.

8. Лисейкин В. Д., Шокин Ю. И., Васева И. А., Лиханова Ю. В. Технология построения разностных сеток. Новосибирск: Наука, 2009.

9. Liseikin V. D. A computational differential geometry approach to grid generation. Berlin: Springer, 2004.

10. Liseikin V. D. A computational differential geometry approach to grid generation. 2nd ed. Berlin: Springer, 2007.

11. Liseikin V. D. Grid generation method. 2nd ed. Berlin: Springer, 2010.

12. Vaseva I. A., Liseikin V. D., Likhanova Yu. V., Morokov Yu. N. An elliptic method for construction of adaptive spatial grids // Russ. J. Numer. Anal. Math. Modelling. 2009. Vol. 1. Р. 65-78.

13. Daniel Appelo, N. Anders Petersson. A Stable Finite Difference Method for the Elastic Wave Equation on Complex Geometries with Free Surfaces // Communications in Computational Physics. 2009. Vol. 5. No. 1. Р. 84-107.

14. Mikhailenko B. G. Spectral Laguerre method for the approximate solution of time dependent problem // Appl. Math. Lett. 1999. Vol. 12. P. 105-110.

15. Демидов Г. В., Мартынов В. Н. Пошаговый метод решения эволюционных задач с использованием функций Лагерра // Сибирский журнал вычислительной математики. 2010. Т. 13, № 4. C. 413-422.

16. Фадеев Д. К., Фадеева В. Н. Вычислительные методы линейной алгебры. M.: Физмат-гиз, 1960.

17. Reinders James, Jeffers James. Intel Xeon Phi Coprocessor High Performance Programming. Morgan Kaufman Publ. Inc, 2013.

18. KlemmMichael. Programming for the Intel Xeon Phi Coprocessor. Intel Corporation, 2013.

19. Rezaur Rahman. Intel Xeon Phi Coprocessor Architecture and Tools: The Guide for Application Developers. Apress, 2013.

20. Parallel Programming and Optimization with Intel Xeon Phi Coprocessors. Colfax, 2013

21. McCoolMichael, Reinders James, Robison Arch. Structured Parallel Programming: Patterns for Efficient Computation. Morgan Kaufman Publ. Inc, 2013.

Материал поступил в редколлегию 16.10.2018

P. A. Titov

Institute of Computational Mathematics and Mathematical Geophysics SB RAS 6Academician Lavrentiev Ave., Novosibirsk, 630090, Russian Federation

[email protected]

ELASTIC WAVES MODELING IN MEDIA WITH COMPLEX FREE SURFACE TOPOGRAPHY

Modeling is an integral part of the study of seismic wave propagation processes in various media. The widely used method of dividing a region into rectangular cells has a drawback: with a complex free surface structure (for example, a mountain), effects arise when a wave is reflected from this surface, due to the fact that the boundary is approximated by a step function for the numerical solving of the problem. In this paper, we propose a different approach to the discretization of the domain: the construction of a curvilinear grid that agrees well with the geometry of the free surface. An algorithm for numerical simulation based on a numerical solution of the linear 2D-system of the elasticity, written in displacements, using a curvilinear grid and the stepwise Laguerre method, is proposed. Results of the simulation are presented. Two parallel versions of the algorithm were also implemented, calculations were carried out on different multicore systems (on the classical multiprocessor architecture, as well as on the architecture using the Intel Xeon Phi co-processors - "RSC PetaStream"). Comparison tests of accelerations on different architectures are presented. Also features of the algorithm parallelization are described.

Keywords: 2D, elasticity, elastic waves, curvilinear surface, curvilinear grids, Laguerre method, Intel Xeon Phi, PetaStream, MPI, OpenMP.

References

1. Liseikin V. D., Rychkov A. D., Kofanov A. V. Adaptive grid technology for applications. Novosibirsk, Novosibirsk State University, 2011, 160 p. ISBN 978-5-94356-981-4 (in Russ.)

2. Khakimzyanov G. S., Shokin Yu. I. Lectures on finite difference schemes on moving grids. Problems for partial differential equations with two space variables. Kazan, Publishing center of KazNU named after Al Farabi, 2006. (in Russ.)

3. Vaseva I. A., Kofanov A. V., Liseikin V. D., Yu. V., Kharitonchik A. M. Constructive technology of spatial adaptive grids. Journal of Computational Mathematics and Mathematical Physics, 2010, vol. 50, № 1, p. 1-19. (in Russ.)

4. Vaseva I. A., Liseikin V. D. Finite element method for adaptive grids construction Computational Technologies, 2011, vol. 16, № 6, p. 30-35. (in Russ.)

166

n. A. Thtob

5. Glasser A., Liseikin V. D., Kitaeva I. A. Differential grids properties controlling via monitor metric. Journal of Computational Mathematics and Mathematical Physics, 2005, vol. 45, № 8, p. 1466-1483. (in Russ.)

6. Glasser A., Liseikin V. D., Shokin Yu. I., Vaseva I. A., Likhanova Yu. V. Differential grids construction via Beltrami and diffusion equations. Novosibirsk, Nauka, 2006. (in Russ.)

7. Liseikin V. D., Morokov Yu. I., Vaseva I. A., Likhanova Yu. V. Some aspects of adaptive grids construction. Journal of Computational Mathematics and Mathematical Physics, 2008, vol. 48, № 9, p. 1638-1658. (in Russ.)

8. Lisekin V. D., Shokin Yu. I., Vaseva I. A., Likhanova Yu. V. Differential grids construction technology. Novosibirsk, Nauka, 2009. (in Russ.)

9. Liseikin V. D. A computational differential geometry approach to grid generation. Berlin, Springer, 2004.

10. Liseikin V. D. A computational differential geometry approach to grid generation. 2nd ed. Berlin, Springer, 2007.

11. Liseikin V. D. Grid generation method. 2nd ed. Berlin: Springer, 2010.

12. Vaseva I. A., Liseikin V. D., Likhanova Yu. V., Morokov Yu. N. An elliptic method for construction of adaptive spatial grids. Russ. J. Numer. Anal. Math. Modelling, 2009, vol. 1, p. 6578.

13. Daniel Appelo, N. Anders Petersson. A Stable Finite Difference Method for the Elastic Wave Equation on Complex Geometries with Free Surfaces. Communications in Computational Physics, 2009, vol. 5, no. 1, p. 84-107.

14. Mikhailenko B. G. Spectral Laguerre method for the approximate solution of time dependent problem. Appl. Math. Lett, 1999, vol. 12, p. 105-110.

15. Demidov G. V., Martynov V. N. Stepwise method for solving evolution problems using Laguerre functions. Siberian Journal of Numerical Mathematics, 2010, vol. 13, № 4, p. 413-422. (in Russ.)

16. Fadeev D. K., Fadeeva V. N. Computational methods of linear algebra. Moscow, Fizmatgiz, 1960, p. 214-220. (in Russ.)

17. James Reinders, James Jeffers. Intel Xeon Phi Coprocessor High Performance Programming. Morgan Kaufman Publ. Inc, 2013.

18. Michael Klemm. Programming for the Intel Xeon Phi Coprocessor. Intel Corporation, 2013.

19. Rezaur Rahman. Intel Xeon Phi Coprocessor Architecture and Tools: The Guide for Application Developers. Apress, 2013.

20. Parallel Programming and Optimization with Intel Xeon Phi Coprocessors. Colfax, 2013.

21. Michael McCool, James Reinders, Arch Robison. Structured Parallel Programming: Patterns for Efficient Computation. Morgan Kaufman Publ. Inc, 2013.

Received 16.10.2018

For citation:

Titov P. A. Elastic Waves Modeling in Media with Complex Free Surface Topography. Vestnik NSU. Series: Information Technologies, 2018, vol. 16, no. 4, p. 153-166. (in Russ.) DOI 10.25205/ 1818-7900-2018-16-4-153-166