Параллельная реализация метода численного моделирования волновых полей в трехмерных моделях неоднородных сред* Текст научной статьи по специальности «Математика»

Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2009, № 6 (1), с. 203-209

УДК 550.344

ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ МЕТОДА ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ВОЛНОВЫХ ПОЛЕЙ В ТРЕХМЕРНЫХ МОДЕЛЯХ НЕОДНОРОДНЫХ СРЕД*

© 2009 г. Д.А. Караваев

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

dmitry1985@ngs.ru

Поступила в редакцию 07.08.2009

Рассматривается численное моделирование распространения упругих волн в трехмерных неоднородных средах, характерных для грязевулканических структур. Численное моделирование проводится на основе конечно-разностного метода с применением метода поглощающих границ. Предлагаются способы распараллеливания программы для численных расчетов средствами языков параллельного программирования MPI и OpenMP. Представлены результаты численного моделирования в виде сейсмотрасс и мгновенных снимков волнового поля.

Ключевые слова: параллельная программа, упругие волны, грязевые вулканы, 3D численное моделирование, неоднородная среда.

Введение

В настоящее время актуальной проблемой является задача генезиса грязевых вулканов, поскольку у исследователей нет единого мнения о происхождении грязевулканических образований. Грязевые вулканы Керченско-Таманской области заслуживают внимательного научного изучения, так как их корни уходят на глубины до 8-10 км и дают интереснейшую информацию для геологов. Нужно отметить, что некоторые грязевые вулканы в Таманской грязевулканической провинции находятся в непосредственной близости от населенных пунктов и представляют реальную опасность для населения в моменты интенсивных извержений. Поэтому создание адекватных моделей происходящих в вулканах этого типа процессов является важной и актуальной задачей.

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

Результаты первых вибросейсмических исследований грязевого вулкана Шуго [1] показа-

* Статья рекомендована к печати программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2009» (http://agora.guru.ru/pavt).

ли, что волновое поле имеет сложное строение и для интерпретации полученной информации и проведения новых полевых экспериментов целесообразно провести математическое моделирование. Проведение численных экспериментов по распространению упругих волн в трехмерных моделях сред, «близких» к исследуемому объекту, позволит выбрать приемлемую расстановку системы возбуждения и регистрации при проведении полевого эксперимента, а также поможет при интерпретации сейсмической информации.

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

Разработанная параллельная программа предназначена для численного моделирования распространения упругих волн в трехмернонеоднородных упругих средах с использованием конечно-разностного метода.

Постановка задачи

Численное моделирование распространения сейсмических волн в сложно построенных уп-

ругих неоднородных средах проводится на основе полной системы уравнений теории упругости с соответствующими начальными и граничными условиями. Данная постановка задачи представлена в терминах вектора скоростей (и, V, Ж) смещений и тензора напряжений (а , а , а , а , а , а ):

V ХХ ’ уу? 22 ’ ХУ хг У2

рЖ = 8^ ^ +80^ + ^ ц, х, у, 2)>

д дх ду дг

д¥ даху да ^ дау2

Р----=-------- +-(?, X, У, 2),

д? дх ду дг

дЖ да„ да да_ „ , .

Р----= —— + —— + —— + К (?, X, У, г),

д? дх ду дг

дахх ди .дУ .дЖ

—— = (X + 2ц)----------------+ X— + X-,

дt дх ду дг

да}у ди . ,д7 ,дЖ

—— = X--------+ (X + 2ц)-----+ X----,

д? дх ду дг

да ди дУ дЖ

да^ = Х— + Х— + (Х + 2ц)—, дt дх ду дг

да хг ди дЖ л

—— = И(------+------),

д/ дz дх

да д¥ дЖ,

— = й(— +--------),

д? дг ду

даху ди дГл

------ = К— + —)

д? дУ дх

(X, ц - параметры Ламе) с начальными

и(х, у, г, ?) |,=0 = 0, V(х, у, г, ?) |(=0 = 0 ,

Ж(х, у, г, t )|(=0 = О и граничными условиями:

а | п = 0, а | п = 0 , а | п = 0.

хг I г=0 ’ У? 1?=0 ’ гг 12=0

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

^ ^ ^ ^

р(I, х, у, г) = р г + р у + р к.

Например, для источника типа «вертикальная сила» получим следующее представление:

^(I, х, у, г) = 5(х - х0 )5(у - у0 )5(г - z0)/() к, где (х0, у0, z0) - координаты источника.

тов (X, Х + 2ц, ц, р, которые могут иметь разрывы), участвующих в разностной схеме, проводится на основе интегральных законов сохранения [6]. Для этого обращается матрица Гука. По формулам, приведенным в работе [6], происходит пересчет новых значений (аппроксимирующих интегральный закон сохранения), затем матрица снова обращается, и получаем матрицу с коэффициентами, которые используются для разрешения схемы.

Конечно-разностная схема имеет второй порядок аппроксимации по времени и пространству [6]. Общий вид конечно-разностной схемы следующий:

П+1 П 1 1

и 1 и 1 п+— п+

Р ук +Р г-1 ук г 2 У * 2 У (О 2 ххух О 2 ххг-1 ) +

2 т Ах

(гъ 2 .1.1 _ ^ 2 .1.14

(О хуг--у+-х О хУг~—У~~х )

Ау

1 1

п+— п+----

2.1.1 «.2.1.1'

(г: 2 . 1 . 1_гт 2 . 1 . 14

Х21 1Х +- Х21-IX--)

, Г 2 2 2, -гп

+ + J хух ,

&

V п+\ - Vп 1

рт +ра-1к ч-^ ч--ъ

т

11 +_ и+_

2 1 1 _______ 2 11

хул !-----X ^ XVI-!-X

2 2 2 2

Ах

11

п+— п+

+ (ст УУУх — СТ хуу-1х ) +

Ау

11

п+— п+

(г: 2 ..1 1 _2 ..1 14

(0 у21] х+- О У21]--х--)

2 2 п

А

+/п

wn+ll - ™" 1 рцк +рук-1 чк - 2 чк - 2

2

т

11

п+— п+

(гг 2 .1. 1_гг 2 14

Х21Л IX- ^ Х21 IX-)

2 2 2 2

Ах

1 1

п +— п +----

(г? 2 .. 1 1 _ -Р. 2 .. 1 14

У211 +------------------X---- ^ У211----X-----)

I 2 2______________________ 2 2* +

Ау

Метод решения задачи

Метод решения поставленной задачи основан на использовании конечно-разностного метода. Алгоритм построения конечно-разностной схемы предложен в статье [6]. Аналогичный подход был развит в работах отечественных математиков [7]. Расчет сеточных коэффициен-

п+— п +

, (а : " 22ух и 2 хгух-1 ) -п

+ I + •/ 2

А

11

п+— п--

22 Ф ххук ххук

т

= (^ + 2р)т

ук

11 г +— /к г — /к

2 __2__

Ах

11

+

+

2

+

+

+

11

+

Здесь Ax, Ay, Az

шаги дискретизации по

V 1 V 1 W 1 W 1

ij +—k ij—к ijk +— ijk—

2 J 2 , ^ 2 2

j Ay + Л j'k Az

1

1

+ (Л + 2ц)

a 2 УУ^ a 2 yy^ _

т

n n

u 1 - u 1

i+— jk i— jk

_ Л 2' 2 ,

_л«- a; +

n n n n

vn 1 - vn 1 wn 1 - wn 1

ij +—k ij —k ijk+— ijk—

2___________ 2 + л 2 2

ijk д + Л ijk

Ay

Az

11

n+— n-

22 a zzijk — a zzijk

т

un 1 — un 1 Vn 1 — Vn 1

i+— jk i— jk ij+—k ij—k

j+ л,„ j2 j2

ijk

Ax

+ (Л + 2ц)

ijk

Ay

-/ 1 W-7 1

¡!k +— г/k—

22

ijk

Az

11 n+— n—

2 1 1 2 1 1 a xzi— jk — a xzi— jk—

2 2 2 2

т

u 1 u 1 w 1 w 1 i— jk i— jk—1 ijk — i—1 jk—

_ c44 ( 2 a 2 +-----2-),

i—j jk—j Az Ax

11

n+— n

a 2 .11 —a 2 .11

a yzi —k— a yzi-k—

2 2________ 2 2 _

т

n n n n

V 1 — V 1 w 1 — W 1

ij —k ij—k—1 ijk — ij—1k—

_ c55 1 1 ( 2 , 2 +---------4-----),

ij—2k—2 Az

11 n+— n—

rr 2 . 1 . 1 , _ rr 2 .1.1,

a xyi— j —k a xyi— j —k

2 2 2 2

Ay

_ c66 1.1 (■

* 1 u 1 V 1 V 1

i — jk i----j—1k ij—k i—1 j—k ґл\

2 ____________2 + 2__________________________ 2 ) (1)

i—2 j —jk ' Ay

Ax

Пример расчета взвешенного коэффициента с66 для расчета а :

,1,11 1 1 ЛЛ-1

сбб 1 1 - (—( I I I )) .

!-2 2 * И ¡к Н-;-1 ¡к Ну-1 к Иг-1 ¡-1к

По аналогичным формулам определяются остальные коэффициенты с55 и с44.

Критерий устойчивости данной схемы [6]:

1

т <-

Vp

i m

1 1 1

ІШКЛІ Ax2 Ay2 Az2

пространственным переменным, т - шаг дискретизации по времени; Vpmax - максимальная скорость распространения упругих волн.

Общая схема вычислений выглядит следующим образом. Вначале, на первом полушаге по времени определяются компоненты вектора скорости смещения, затем на втором полушаге по времени по формулам находятся нужные компоненты напряжений. Далее насчитываются новые компоненты образов вектора скорости смещения на новом полушаге по времени.

Поглощающие границы

В связи с тем, что область расчета ограничена, необходимо использовать поглощающие границы. В полной мере данный подход (Perfectly Matched Layers) изложен в работе [8]. В поглощающей зоне происходит пересчет основных искомых величин по конечно-разностным формулам (2), а вне зоны поглощения расчет ведется по формулам (1).

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

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

x y z

а = а + а + а , а =

xx xx xx xx yy

x y z x y z

= а +а y +а ,а =а +а y +а .

УУ yy yy ’ zz zz zz zz *

x y x z

a = a + a , a = a + a , a =

xy xy xy xz xz xz yz

yz

= a y +a .

yz yz *

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

Тогда применение метода [8] дает уравнения для расчета в зоне поглощения:

G xxijk G xxijk j % G xxijk G xxijk

1Л і 1Л і

i-i—jk i—jk

= (X + 2^)jk—-— ----------------- -----+ f "xj,

J Ax

11------------------11 ,,n+— vn- vn+— v"

_ y 2 z^y 2 z^y 2 r^y 2

O xxijk O xxijk + ^ y O xxijk O’ xxijk

T

nn v-1 k - v" 1 k

i +— k и— k

22

- K ijk

Ay

+

т

1 1 — ,.n--------

11

n+— ,.n-

2

T

2

11-------------------------------11 zn+----------------------------- z n- z n+ z n--

C 2 xxijk C 2 xxijk + ^ z C 2 xxijk C 2 xxijk

T

2

(2)

УУ 1 УУ 1

ijk +— ijk—

_ 4 2 2

_ Л ijk Az

В этих уравнениях величины йх, йу, й2 есть коэффициенты поглощения вдоль соответствующих направлений пространственных переменных. Определив из уравнений значения компо-

x у z

нент axx , axx , Gxx

значение компоненты

схх получаем суммированием всех найденных

х , У , 2 “

компонент ахх = ахх + ахх + ахх для каждой точки конечно-разностной схемы в зоне РМЬ. Для оставшихся конечно-разностных уравнений проделываются аналогичные действия для получения значений искомых величин в зонах поглощения.

Параллельная реализация программы численного моделирования

В данной работе решается трехмерная динамическая задача моделирования распространения упругих волн. Как уже было замечено, эта задача является ресурсоемкой, поэтому логично проводить столь массивные вычисления на современных параллельных вычислительных комплексах. Необходимо отметить, что данная параллельная реализация не предполагает использования различных типов ускорителей или сопроцессоров, обладающих нестандартной архитектурой. Для создания параллельного варианта программы рассматривалось несколько подходов. Первый подход к распараллеливанию алгоритма заключается в том, что исходная область расчетов дробится на кубы в количестве, равном количеству вычислительных узлов, имеющемуся в распоряжении пользователя; второй основан на разбиении области на слои вдоль координаты Oz. В первой ситуации резко возрастает количество обменов информацией между гранями, ребрами и вершинами соприкасающихся кубов, однако объем информации, которым обмениваются «соседи», меньше, чем объем информации, который передается соседним узлам при разбиении расчетной области на слои. В данном варианте программы осуществлен второй подход, в первую очередь из-за простоты реализации. В этом случае количество слоев определяется количеством свободных вычислительных узлов. На основе этой реализации созданы две параллельные программы: в первой для распараллеливания используется

только OpenMP, во второй используется комбинация возможностей MPI и OpenMP. В таком случае предлагается проводить обмен информацией между соседними слоями через MPI, а внутри каждого слоя проводить параллельные вычисления используя OpenMP. Приведенные в данной работе расчеты проводились в Институте вычислительной математики и математической геофизики на кластере НКС-30Т, использовалось 2 вычислительных узла по 8 ядер на каждом. Сравнение времени работы этих двух параллельных программ не проводилось, хотя это является интересным и необходимым моментом, и в ближайшее время эти исследования будут выполнены.

Реализация «построителя» трехмерной модели упругой среды

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

В нашем случае предполагается, что задана крупноблочная модель среды, составленная из параллелепипедов, в вершинах которых задаются параметры среды ( Vp,Vs, р ). Физический смысл данных параметров - скорости продольных и сдвиговых волн. Эти параметры являются непрерывными внутри каждого блока. Разрывы проходят только по граням соседних параллелепипедов. Далее методом конечных элементов происходит интерполяция параметров среды на более «мелкую» расчетную сетку.

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

Разработанный построитель модели позволяет конструировать сложные модели сред, близкие к реальным объектам исследования.

Результаты численного эксперимента

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

Для исследования структуры волнового поля, получаемой при проведении полевых экспериментов на грязевых вулканах в Керченско-Таманской области, были проведены тестовые

расчеты для различных моделей сред, одна из которых представлена далее.

Модель упругой среды

Моделировалась среда, состоящая из трех слоев, два из которых имеют куполообразную форму, и цилиндрической подобласти (рис. 1). Размеры среды в направлении координат X - 3.8 км, У - 2.4 км, Z - 2.2 км. Параметры:

- верхний слой - ¥р = 1.0 км/с, ¥8 = 0.5 км/с, р = 2 г/см3;

- средний слой - ¥р = 1.5 км/с, ¥ї = 0.7 км/с, р = 2 г/см3;

- нижний слой - ¥р = 2.0 км/с, ¥ї = 1.0 км/с, р = 2 г/см3;

- цилиндр - ¥р = 1.5 км/с, ¥ї = 0.4 км/с, р = 2 г/см3 , диаметр 0.2 км, координаты оси (2.2 км, 1.2 км);

Система возбуждения

Источник типа «центр давления», несущая частота 6 Гц, близок к свободной поверхности, координаты (0.3 км, 1.2 км, 0.03 км).

Система наблюдения

Вблизи от свободной поверхности расположены 4 линии сейсмоприемников. Схема расположения представлена на рис.1.

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

Заключение

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

С помощью созданной программы моделирования проведены расчеты на многопроцессорных вычислительных системах.

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

Автор выражает благодарность научным руководителям Б.М. Глинскому и В.Н. Мартынову.

Работа выполнена при поддержке: Интеграционный проект СО РАН № 26, Интеграционный проект РАН № 16.5, Программ фундаментальных исследований Президиума РАН № 2.2.

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

1. Глинский Б.М., Собисевич А.Л., Хайретдинов М.С. Опыт вибросейсмического зондирования сложно построенных геологических структур (на примере грязевого вулкана Шуго) // Докл. РАН. 2007. Т. 413, № 3. С. 398-402.

2. Глинский Б.М., Фатьянов А.Г. Численноаналитическое моделирование волновых полей в разномасштабных зонах вулканической деятельности // КВМ-2007. Новосибирск, 18-20 июня 2007.

3. Глинский Б.М. Фатьянов А.Г., Вибросейсмиче-ский мониторинг живущих вулканов // Активный геофизический мониторинг вулканов: Материалы 2-го Межд. симпозиума 12-16 сент. 2005. Новосибирск.

С. 57-61.

4. Глинский Б.М. Фатьянов А.Г. Изучение и мониторинг грязевых вулканов активными сейсмическими методами // Материалы 2-го Межд. симпозиума 12-16 сент. 2005. Новосибирск. С. 52-57.

5. Михайленко Б.Г. Сейсмические поля в сложно построенных средах. Новосибирск, 1988. С. 312.

6. 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). Vol. 2. P. 75-80.

7. Коновалов А.Н. Сопряженно-факторизованные модели в задачах математической физики // Сиб. журн. вычисл. математики / РАН. Сиб. отделение. Новосибирск, 1998. Т. 1. № 1. С. 25-57.

8. Collino F., Tsogka C. Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media // Geophysics. 2001. 66 (1). P. 294-307.

PARALLEL IMPLEMENTATION OF WAVE FIELD NUMERICAL MODELING METHOD IN 3D MODELS OF INHOMOGENEOUS MEDIA

D.A. Karavaev

Numerical modeling of elastic wave propagation in 3D inhomogeneous media typical of mud volcanoes is considered. The numerical modeling is carried out on the basis of the finite-difference method with application of the absorbing boundary method. Parallel languages MPI and OpenMP are proposed to be used for numerical modeling program parallelization. The results of numerical modeling in the form of seismic traces and wave-field snapshots have been presented.

Keywords: parallel program, elastic waves, mud volcanoes, 3D numerical modeling, inhomogeneous medium.