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

УДК 517.97: 539.37

Вычислительный алгоритм для расчета волновых полей в блочных средах на многопроцессорных вычислительных системах

Елена В.Кучунова*

Институт математики, Сибирский федеральный университет, пр. Свободный 79, Красноярск, 660041,

Россия

Владимир М.СадовскиЙ

Институт вычислительного моделирования СО РАН, Академгородок, Красноярск, 660036,

Россия

Получена 11.01.2008, окончательный вариант 15.03.2008, принята к печати 10.04.2008 Разработан параллельный алгоритм для численного решения пространственных задач динамики упругих кусочно-однородных сред блочной структуры с криволинейными поверхностями раздела. Алгоритм основан на методе расщепления по пространственным переменным с применением монотонной БКС-схемы для решения одномерных гиперболических систем. Исследуются вопросы эффективной реализации алгоритма на многопроцессорных вычислительных системах с распределенной памятью. Приводятся результаты апробации комплекса параллельных программ при исследовании модельных задач сейсмики.

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

В данной работе численно решается прямая динамическая задача теории упругости по восстановлению неизвестных полей вектора скорости v(x,t) и симметричного тензора напряжений a(x,t) в некоторой области Q С R3 по заданным начальным данным и внешним воздействиям. Структура среды известна и представлена набором из

m

m разнородных блоков с криволинейными границами: Q = У Qk. Каждый блок Qk ха-

k=l

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

р v = V • а,

(1)

а = A(V • v)I + p(V <g> v + v <g) V).

Здесь р — плотность материала, A и ц — параметры Ламе, определяющие упругие свойства, V — градиент по пространственным координатам, I — единичный тензор. Точка

* e-mail: vek@krasu.ru te-mail: sadov@icm.krasn.ru © Siberian Federal University. All rights reserved

над символом, как обычно, означает производную по времени. Компоненты свертки градиента с тензором напряжений в декартовой системе координат задаются равенством (V • а)і =5^ V і ® — оператор тензорного произведения: (V <8> ч)з = V'.

з

Предполагается, что в 0 задано начальное распределение скоростей и напряжений:

ч(х, 0)= ч0(х), а(х, 0)= а0(х). (2)

Кроме того, известны внешние силы, действующие на части границы 0 на дневной поверхности:

а • п = q, (3)

где п — вектор внешней нормали. На остальной части задаются неотражающие гра-

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

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

Рис. 1. Пример сетки в многоблочной области

Для численного исследования динамических задач механики деформируемых сред, описываемых гиперболическими системами уравнений, широко применяются сеточнохарактеристические методы [1, 2]. В монографии [3] представлено современное состояние этих методов. При построении разностной схемы в задаче (1)-(3) исходную систему

(1) запишем в матричной форме:

д и д и

Ж - ^ Аі дХІ = 0.

і=1

(4)

Здесь и = [V, а]т = («і, «2, «з, ац, а22, азз, а2з, аіз, аі2)т — вектор-функция, со-

ставленная из компонент вектора скорости и тензора напряжений, Аі — матрицы-коэффициенты размерности 9 х 9, произведение которых на вектор и задается равенством

1

л ' ~ єі •а

А- = р

Л(ві • V)/ + ^(ві ® V + V ® Єі)

где ві, в2, вз — единичные взаимно-ортогональные векторы.

Для решения системы (4) применим схему интегрирования на основе метода двуциклического расщепления [4, 5]. При этом расщепление будем производить не в физическом, а в параметрическом пространстве. Исходную систему преобразуем к виду

V

а

ди = "ТХ— А- = "Т ^А-

а* =2^Аід^ Аі = 2^ дхАз.

дх-з=і -1

(5)

С учетом обозначения а® = V • получим следующую формулу для вычисления произведения матриц на вектор:

А-

V

а

1 -— а • а

р

Л(а- • V)/ + м(а- <8> V + V ® а-)

Применяя метод двуциклического расщепления по пространственным переменным к системе (5), приходим к серии из шести одномерных задач на интервале * Є [*п,*п + т]:

ди(1) ~ ди(1)

= Аі-

д*

дЄі ’

ди(2) ~ ди(2)

= ^42-

д*

д^2

ди(з) 4 ди(з)

= Аз-

д*

д£з

ди(4) ~ ди(4)

= Аз-

д*

д£з

ди(5) ~ ди(5)

= А2-

д*

д^2

ди(6) 4 ди(6)

Аі-

д*

дЄі ’

и(1)(*п) = и(*„), * Є (*„,*„ + т/2);

и(2)(*„) = и(1)(*„ + т/2), * Є (*„,*„ + т/2);

и(з)(*„) = и(2)(*„ + т/2), * Є (*„,*„ + т/2);

и(4)(*„ + т/2) = и(з)(*„ + т/2), * Є (*„ + т/2,*„ + т);

и(5)(*„ + т/2) = и(4)(*„ + т), * Є (*„ + т/2,*„ + т);

и(6)(*„ + т/2) = и(5)(*„ + т), * Є (*„ + т/2,+ т).

(6)

Вектор-функция и(6) (£„ + т) представляет собой искомое решение на временном слое *п + т.

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

І = (7)

где Q - матрица, составленная из линейно-независимых собственных векторов, Л — диагональная матрица, элементами которой служат собственные значения. Для направления £ = £

Л = а1ая{ер|а-|, -Ср|а-|,е8|а-|, с8 |а-|, -с8|а-|, -с8|а-|, 0,0, 0},

где |а-| = -у^(аі)2 + («2)2 + (аз)2, ср = \/(Л + 2й)/р — скорость продольных волн, ся = \/й/р — скорость поперечных волн в среде.

К одномерным системам (7) применяется явная монотонная ЕМО-схема предиктор-корректор с предельной реконструкцией решения (см. [3]). На шаге “предиктор”, исходя из характеристических соображений, по известным значениям на п-м временном слое вычисляют значения решения на границах ячеек. Система (7) после замены неизвестной вектор-функции по формуле w = Qu распадается на девять независимых скалярных уравнений переноса

д w дw

Ж = д£ ,

для решения которых из узла (£-, + т/2) (_?’ = 0,1,..., N) опускаем характеристики

на слой * = (рис. ). Координаты вектора w постоянны вдоль характеристик, поэтому

^П+1/2 (£з) = ^П(£з + Л к і+1/2 т/2),

где Л д,і+і/2 — угловой коэффициент к-й характеристики в ]-й ячейке. Характеристики, проведенные из узла (£-,*п + т/2), могут не попасть точно в расчетные узлы (£-±1/2,*п), в которых определено решение, поэтому применяется процедура предельной реконструкции: сеточная функция с узловыми значениями ^-+1/2 на каждом из отрезков [£-, £-+і] заменяется кусочно-линейным сплайном г4П(£) = аП3'+і/2£ + ^ 1/2,

коэффициенты которого подбирают из условия минимизации суммы модулей скачков сплайна на границах ячеек [3]. В итоге получаем следующий вектор w"+1/2:

+ 1/2 _ I т Пі+1/2 к32+1/2 (£І + 1 - ) + 2 Л к і'+1/2«Пі+і/2, если Л к і+1/2 > 0,

д з — 1 ап -

1тП3'-1/2 + 12 (£з — £з-і) + 2 Л к,І-1/2 - — і/2, если Л к і — і/2 < °.

После этого восстанавливается само решение:

и”+1/2 = Q-1w”+1/2. (8)

Формулу (8) нельзя использовать в приграничных ячейках, поскольку в нее входят величины т П^±1/2 в точках, которые находятся за пределами области интегрирования. Корректное решение задачи на границе требует постановки граничных условий в количестве, равном числу уходящих характеристик [2]. Для левой границы области

Рис. 2. Сеточно-характеристический метод расчета на шаге “предиктор”

уходящими являются характеристики, отвечающие отрицательным собственным значениям, для границы — положительным (рис. 2). В качестве недостающих условий на дневной поверхности области Q выступают граничные условия (3). На остальной части поверхности ставим искусственные граничные условия, основываясь на принципе неот-ражения волн. Следуя этому принципу, полагаются равными нулю граничные значения характеристических переменных на уходящих характеристиках, что эквивалентно отсутствию отраженных волн в одномерной задаче.

Для сохранения свойств консервативности схемы на криволинейных сетках, когда матрицы-коэффициенты Aj зависят от координат, аппроксимация одномерных систем уравнений на шаге “корректор” производится с помощью интегро-интерполяционного метода. В результате интегрирования системы (5) по криволинейной ячейке в физической области с последующим применением формулы Грина получается равенство

о— 1 6

Ж = 7к(пк^ + nk A2 + nk Аз) —(9)

fc=i

в котором nj — направляющие косинусы внешней нормали, w — объем ячейки, и — среднее интегральное значение вектора-решения, а индекс к относится к граням ячейки, в частности Yk — площадь соответствующей грани. Далее сумма в правой части (9) разбивается на три пары слагаемых по противоположным граням, каждая из которых отвечает аппроксимации производных по пространственным переменным в одномерных системах (6). Например, слагаемые, относящиеся к координатным поверхностям £2,£з,

~ ди

дают аппроксимацию выражения Ai ——.

d£i

Заметим, что применяемая схема решения одномерных задач является вариантом сеточно-характеристической схемы второго порядка точности. Метод двуциклического расщепления сохраняет второй порядок. Кроме того, двойной пересчет одной и той же системы уравнений на третьем и четвертом этапах расщепления гарантирует устойчивость метода при выполнении одномерного условия Куранта-Фридрихса-Леви.

Изложенный выше численный алгоритм реализован комплексом программ на алгоритмическом языке Fortran-90 с использованием библиотеки обмена сообщениями MPI (Massage Passing Interface). Подробно принципы параллельной реализации алгоритма изложены в [8]. Технология распараллеливания изложенного выше вычислительного алгоритма основана на разделении исходной области по узлам в вычислительной системе, исходя из требования равномерной загрузки узлов.

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

Р(і) = Ро бій (ші), если і < Т, (10)

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

Задача 1. Пусть изотропная упругая среда занимает прямоугольную область О = {х Є Д3| 0 < жі < Н, — Н < Х2 < Н, — Н < жз < Н} в декартовой системе координат. Параметры среды: р = 7850 кг/м , ср = 6000, е8 = 3210 м/с. В точке хо = (0, 0, 0) приложена точечная импульсная нагрузка (10), действующая перпендикулярно к плоскости ж і = 0. Размеры расчетной области: 200 х 400 х 400 м3. Расчеты проводили на сетке из 200 х 400 х 400 узлов. Шаг по времени выбирали с учетом условия Куранта-Фридрихса-Леви, он составил 0.0016 с. Время действия нагрузки Т = 0.005 с, то есть около четырех шагов по времени. Величина амплитуды колебаний Ро = 250 МПа. Среднее время расчета одного шага по времени на 32-х вычислительных узлах кластера МВС-1000 — 2 минуты.

Рис. 3. Сейсмограмма вертикальной компоненты вектора перемещения и і

Рассмотрим полученную волновую картину в вертикальном сечении области, проходящем через точку приложения нагрузки хо. На рис. 3 приведена сейсмограмма вертикальной компоненты вектора перемещения иі. На рис. 4 представлена визуализация и і в фиксированный момент времени 0.016 с. Здесь введены следующие цифровые обозначения: 1 — продольная волна, 2 — поперечная волна, 3 — головная поперечная волна,

4 — поверхностная волна Релея.

Результаты соответствуют точному решению задачи для случая плоского деформированного состояния: видны продольная и поперечная волны, менее интенсивные головные волны и волны Релея, которым на рисунках соответствуют движущиеся точки на поверхности ж і =0.

Задача 2. Рассмотрим модель двухслойной среды с однородным изотропным верхним слоем и скачком акустической жесткости на границе с подстилающим полупространством. Под однородной покрывающей средой с параметрами р=500 кг/м3, ср=500

( !!!"""$ !!!"""' !!!"?"$""# !!!"""& !!!"""( !( !!!

Рис. 4. Визуализация вертикальной компоненты вектора перемещения п\ в вертикальном разрезе в момент времени t=16 мс

и cs=300 м/с расположена вторая среда с параметрами ро=750 кг/м3, cpj=750 и е°=450 м/с (рис. 5). Горизонт, разделяющий две среды, представляет собой плоскость, наклоненную по линии AC под углом у>1 =5°, а по линии DB — под углом ^2 = 7°. Заметим, что на границе раздела сред выполняются условия pcp = pocpj и pcs = роС?, гарантирующие образование отраженных волн.

Рис. 5. Модель области с наклонным горизонтом раздела двух сред

Размеры расчетной области: АВ = АВ = 600, АА2 = 400 м. Глубина залегания наклонного горизонта: АА1 = 175, ВВ1 = 87, 5, СС1 = 100 , ВВ1 = 187, 5 м. Внешний импульс на среду действует в точке О. Глубина залегания границы раздела в точке О составляет ОО1 = Н = 156, 25 м. Величина приложенной нагрузки — Ро=250 МПа, время действия — 0,01 с, радиус действия — 3.5 м. Число узлов сетки — 200 х 200 х 200, временной шаг — т = 0.002 с. Число используемых процессоров — 8=2х2х2, количество

расчетных ячеек, приходящихся на один процессор, — 106. Число шагов по времени — 1000, среднее время расчета одного шага — 2.63 минуты.

На рис. 6 представлены сейсмограммы вертикальной составляющей вектора перемещения вдоль профилей ОС и ОБ. Здесь введены следующие обозначения: 1 — падающая Р-волна; 2 — падающая 8-волна; 3 — поверхностная волна Релея; 4 — отраженная РР-волна; 5 — отраженная Р8-волна; 6 — преломленная РР-волна; 7 — преломленная Р8-волна; 8 — головная волна, образованная РР-волной; 9 — вторичные отраженные волны от верхней границы области; 10 — отраженная 88-волна; 11 — отраженная 8Р-волна; 12 — преломленная 8Р-волна; 13 — преломленная 88-волна; 14 — головная волна, образованная преломленной 8Р-волной; 15 — головная волна, образованная преломленной 88-волной.

Рис. 6. Сейсмограмма вертикальной составляющей вектора перемещения вдоль трассы OC (а) и вдоль трассы DB (б) для модельной задачи 2

При прохождении прямых P- и S-волн через границу раздела сред образуются преломленные и отраженные волны. Годограф отраженной PP-волны вдоль профиля AC характеризуется уравнением: t(r) = — \Jr‘2 — 2Hr sin 2y>i + 4H2 cos2 и

Cp

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

-2ІТ-

Рис. 7. Фрагменты визуализации вертикальной компоненты вектора перемещения в моменты времени 700мс (а) и 900мс (б)

подъема горизонта на расстояние ri = H sin2wi=27 м, время прихода отраженной волны ti = 2Hcos wi/cp = 0.623 с. Эта точка отмечена на сейсмограммах Pi. Вдоль профиля DB уравнение годографа отраженной волны усложняется тем, что угол горизонта зависит от r, расстояния до центра сейсмограммы:

tan wi + (r/l) tan ^ 2

w = arctan-------, =-----. (11)

^1+ r2/l2

Годограф отраженной волны вдоль профиля DB:

t(r) = — уr2 + l2 — 2Нл/12 + r2 sin 2w + 4H2 cos2 w .

Cp

Годограф отраженной SS-волны вдоль линии AC описывается уравнением

t(r) = —\Jr2 — 2Hr sin 2w + 4H2 cos2 w .

Cs

Годограф отраженной SS-волны вдоль профиля DB:

t(r) = —у r2 + l2 — 2Hv^2 + r2 sin2w + 4H2 cos2 w .

Cs

При переходе продольной P-волны через границу раздела двух сред кроме отраженной продольной волны образуется вторая отраженная волна — поперечная PS-волна. Аналогично при прохождении поперечной S-волны кроме отраженной поперечной SS-волны образуется отраженная продольная SP-волна. Эти волны обозначены на сейсмограммах цифрами 5 и 11 соответственно.

Критический угол i, при падении под которым луч продольной P-волны преломляется под углом в 90° к нормали, определяется по закону Снеллиуса: sini = cp/cpp. При

- 21В -

этом преломленная продольная волна начинает распространяться вдоль границы раздела сред. Фронт проходящей продольной PP-волны, скользя вдоль границы раздела, возбуждает в верхнем слое колебания, которые вызывают появление головной преломленной волны. Одной стороной фронт головной волны касается фронта отраженной из критической точки волны, другой примыкает к фронту скользящей преломленной волны. На рис. 7 представлены фрагменты визуализации вертикальной компоненты вектора перемещения для случаев падения продольной P-волны под углом, меньшим критического угла i (а) и большим критического (б). При этом на рис. 7 (б) прослеживается образование головной продольной волны. Ее годограф для профиля AC имеет вид

t(r) = —(r sin(i ^ wi) + 2Hcos wi cosi), (12)

Cp

где знак “ —” берут для поднимающегося горизонта, а знак “ +” — для опускающегося. Точка выхода головной волны на поверхность (точка P2 на сейсмограмме) расположена от источника на расстоянии: r2 = 2Hcos wsin i/cos(i — w), время выхода головной волны на поверхность составляет:

2H cos2 wi cp cos(i — w)

Картина распространения головной поперечной SS-волны аналогична. Уравнение годографа поперечной головной волны имеет вид:

t(r) = —(r sin(i ^ wi) + 2Hcos wi cosi).

cs

Точка выхода (точка P3 на сейсмограмме) на поверхность расположена от источника на расстоянии гз = 2H cos wi sin i/cos(i — wi) = 261 м, время выхода составляет t3 = 2H cos2 wi/(cs cos(i — w)) = 1.475 c. Головная поперечная волна выходит на поверхность в тех же точках, что и головная продольная волна, отличается только время. Таким образом, результаты расчетов полностью соответствуют волновой картине, полученной по закону геометрической оптики.

В работе исследуется вопрос эффективности распараллеливания алгоритма. Понятие эффективности рассматриваем как отношение ускорения к количеству процессоров, на котором оно достигнуто [6, 7], в процентах. Анализ эффективности распараллеливания на модельных задачах показывает 80-90 %, что служит довольно хорошим результатом, с учетом неизбежных затрат времени на необходимые пересылки и запись результатов в двоичные файлы.

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

Работа выполнена при финансовой поддержке Красноярского краевого фонда науки (шифр гранта 17G036), РФФИ (код проекта 08-01-00148) и Комплексной программы

фундаментальных исследований Президиума РАН № 14 “Фундаментальные проблемы информатики и информационных технологий”.

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

[1] И.Б.Петров, А.С.Холодов, Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом, Ж. вычисл. матем. и матем. физ., 24(1984), №5, 722-739.

[2] К.М.Магометов, А.С.Холодов, Сеточно-характеристические численные методы, М., Наука, 1988.

[3] А.Г.Куликовский, Н.В.Погорелов, А.Ю.Семенов, Математические вопросы численного решения гиперболических систем уравнений. М., Физматлит, 2001.

[4] Г.И.Марчук, Методы расщепления. М., Наука, 1988.

[5] О.В.Садовская, Метод сквозного счета для исследования упругопластических волн в сыпучей среде, Ж. вычисл. математики и матем. физики, 44(2004), №10б, 19091920.

[6] В.В.Воеводин, Вл.В.Воеводин, Параллельные вычисления. Спб., БХВ-Петербург, 2002.

[7] Дж.Ортега, Введение в параллельные и векторные методы решения линейных систем. М., Мир, 1991.

[8] Е.В.Кучунова, О.В.Садовская, В.М.Садовский, Комплекс прикладных программ для численного решения пространственных задач динамической теории упругости на многопроцессорных вычислительных системах, Избранные материалы Четвертой школы-семинара ’’Распределенные и кластерные вычисления”, Красноярск, ИВМ СО РАН, 2005, 159-172.

Computational Algorithm for Calculation of Wave Fields in Block Media on Multiprocessor Computers

Elena V.Kuchunova Vladimir M.Sadovskii

We develop a parallel algorithm for finding numerical solutions of 3D dynamic problems of the theory of elasticity in piecewise homogeneous block media with curvilinear interfaces. This algorithm is based on the procedure of splitting with respect to spatial variables, in each stage of which the monotone ENO-scheme is used to find solutions of 1D hyperbolic systems. The questions of effective realization of the algorithm on multiprocessor computers with distributed memory are considered. The developed package of parallel programs is used for investigation of model problems in seismics. Keywords: dynamics, elasticity, seismic waves, finite-difference method, decomposition, parallel algorithms.

- 22G -