CC BY
59
УДК 550.834
М.А. Белоносов
ИНГГ СО РАН, Новосибирск
ПРИМЕНЕНИЕ ИНТЕГРАЛЬНОГО ПРЕОБРАЗОВАНИЯ ЛАГЕРРА И МЕТОДА ДЕКОМПОЗИЦИИ ОБЛАСТИ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ СЕЙСМИЧЕСКИХ ВОЛНОВЫХ ПОЛЕЙ
В настоящей работе для численного моделирования процессов формирования и распространения волн в упругих неоднородных средах используется подход, основанный на применении интегрального преобразования Лагерра и последующей декомпозиции области. Несомненным преимуществом этого преобразования является отрицательная определённость возникающего в результате его применения дифференциального оператора, что существенно облегчает использование в дальнейшем метода альтернирования по Шварцу на основе декомпозиции расчётной области. Декомпозиция осуществляется таким образом, чтобы появляющиеся при этом элементарные подобласти могли бы быть загружены на отдельный процессорный элемент. Возникающие при этом системы линейных алгебраических уравнений решаются на основе применения LU разложения соответствующих матриц, которое выполняется для каждой из подобластей только один раз.
Данный подход реализован для вычислительных систем с параллельной архитектурой. Приводятся результаты численных расчётов для реалистичной модели одного из районов Северного моря. Исследована зависимость числа итераций метода альтернирования по Шварцу в зависимости от количества подобластей.
M.A. Belonosov
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS (IPGG)
Acad. Koptyug av. 3, Novosibirsk, 630090, Russian Federation
IMPLEMENTATION OF THE INTEGRAL LAGUERRE TRANSFORM AND DOMAIN DECOMPOSITION METHOD FOR NUMERICAL SIMULATION OF ELASTIC WAVE FIELDS
An approach based on application of the integral Laguerre transform and subsequent domain decomposition for simulation of particle formation and propagation of elastic waves in heterogeneous mediums is proposed. The apparent advantage of this transformation is that we deal with negative definite differential operator, which is much easier to use in future Schwarz alternative method based on decomposition of the computational domain. Decomposition is carried out in such a way that an appearing elementary subdomain could be loaded on a separate processing element. Systems of linear algebraic equations, emerging this way, are solved by application of LU factorization of the matrix, which is performed for each subdomain, only once.
This approach is implemented for computer systems with parallel architecture. The results of numerical calculations for a realistic model of one region of the North Sea is presented. The dependence of the number of iterations of Schwarz alternative method, depending on the number of subdomains is studied.
Введение
Численное моделирование сейсмических волновых полей для реалистичных трёхмерных моделей невозможно без привлечения современных вычислительных систем с параллельной архитектурой.
Наиболее распространённым способом организации таких вычислений является декомпозиция расчётной области на элементарные подобласти с последующим использованием явных конечно-разностных схем. Обладая неоспоримыми преимуществами, такими как простота и эффективность реализации, этот подход имеет и ряд недостатков, к которым относится необходимость выполнения обменов между соседними процессорными элементами на каждом шаге по времени, выбор шага дискретизации по времени исходя из максимальной скорости распространения волн, а шага дискретизации по пространству, наоборот, исходя из минимального её значения. Тем самым моделирование волновых полей в средах с резкими контрастами становится весьма ресурсоёмким.
Необходимость повышения уровня извлечения углеводородов из сложноустроенных трещиновато-кавернозных карбонатных коллекторов Восточной Сибири привела к постановке нового класса задач -моделирования процессов возникновения и распространения сейсмических волн, рассеянных на микронеоднородностях в таких коллекторах. Естественно, что для их описания требуется чрезвычайно мелкая сетка, как по пространству, так и по времени. Использование такой сетки во всём пространстве выдвигает абсолютно нереалистичные требования на вычислительные ресурсы. Альтернативой этому являются сетки с локальным измельчением по времени и пространству. Надо отметить, что реализация такого подхода для явных конечно-разностных схем весьма сложна и связана со значительными накладными расходами ([4], [5], [10]). В то же время, при реализации метода декомпозиции и последующего применения метода альтернирования по Шварцу такое локальное измельчение может быть выполнено именно в отдельно взятой элементарной подобласти, лежащей на специально выделенном для этого процессорном элементе.
Преодоление проблем, описанных выше возможно, если решать задачу с отделением времени. Именно реализации этого подхода посвящена настоящая работа. При этом отделение по времени выполняется посредством интегрального преобразования Лагерра [6] и мы имеем дело с отрицательно определённым дифференциальным оператором. Показывается, что метод альтернирования по Шварцу будет в данном случае весьма эффективен.
Отметим, что использование преобразования Фурье приводит к задаче со знаконеопределённым оператором, что влечёт за собой множество трудностей (см. например [3], [8], [11]).
Преобразование Лагерра и метод Шварца
Рассмотрим интегральное преобразование Лагерра по времени [6]
0
с формулой обращения
да
и(X, 2, ?) = ^ ип (X, г)-(Ы)2 а (Ы), п=0
где а (Ы) - ортонормированные функции Лагерра
1а аы)=
Ип!
а
Ы
(п + а)!
(Ы)2 е 2 Ь<а (Ы).
Здесь п е2 , Ы е Я (Ы >0) , аеЪ (а> 0) , щ (Ы) - классические стандартизованные полиномы Лагерра [9].
Применяя преобразование Лагерра (1) к двумерной системе уравнений теории упругости второго порядка, мы получаем следующую систему уравнений
д_
дх
д_
дх
(Л + 2М)(ып) х +Л(иЦ)
п
д +— дг
и(ип ) х +№(и1 )
п
Ы
рт и1
п
= ( х, г, и°
0
«п-1\ ^ и1 Л
М(и2 )х + Ц(и\ )
п
+
д_
дг
Л(и1) х + (Л + 2ц)(ип)
п
= /? (х, г, и2
Ы_ п
РТ и2
п-К
^ и2 X
(3)
где р - плотность, Ли /и - коэффициенты Ламе, и /2 - правые части. Аналогично тому, как это было сделано в [12], поставим на границе расчётной области условия РМЬ [1].
К полученной системе (3) мы можем применить метод декомпозиции области, основанный на методе альтернирования по Шварцу [2]. Для того чтобы это сделать мы разбиваем расчётную область на несколько пересекающихся подобластей, как показано на рис. 1а. Далее, для уравнений в каждой подобласти применяем конечно-разностную аппроксимацию на сдвинутых сетках и получаем систему линейных алгебраических уравнений. В параллельной реализации число итераций по Шварцу это, вообще говоря, число обменов между смежными подобластями, поэтому их минимизация очень важна. Для системы уравнений (3) число итераций близко к числу подобластей, на которые мы разбиваем расчётную область. Более того, существуют модификации этого метода [7], что число итераций равно числу подобластей.
г
б
Рис. 1. а - декомпозиция области (серым цветом обозначены области, где пересекаются две подобласти, черным - где не пересекаются); б - North Sea
model
Отметим, что, применяя наш подход, мы сталкиваемся с многократным решением систем линейных алгебраических уравнений. Одной из особенностей применения преобразования Лагерра является то, что мы получаем системы с одинаковой матрицей левой части. Благодаря этому мы проделываем LU разложение этой матрицы один раз. Дальнейшее решение системы производится с использованием процедур Math Kernel Library, а именно процедур DSS.
Численный эксперимент
Для демонстрации предложенного алгоритма, основанного на применении интегрального преобразования Лагерра и последующего метода альтернирования по Шварцу, была проведена серия численных экспериментов. В ней использовалась упругая среда, описывающая один из районов Северного моря, представленная на Рис. 1б. Возбуждение волнового поля производилось источником типа «центр объемного расширения» с импульсом Риккера, обладающим доминирующей частотой 30 Гц. Источник располагался в точке с координатами (1600,20). Шаг дискретизации по пространству был взят постоянным и равным 2 м. Для ограничения расчётной области использовался идеально подходящий поглощающий слой (PML - от английского Perefectly Matched Layer). Были взяты следующие значения для параметров функций Лагерра (2):
h = 300, n = 500, а = 5.
Волновое поле рассчитывалось на 3 сек.
Область разбивалась на 12 и 20 подобластей, как показано на Рис. 1а, с шириной перекрытия 60 точек расчётной сетки. Оказалось, что для разбиения на 12 подобластей метод альтернирования по Шварцу сходится за 14 итераций, а на 20 подобластей - за 24. Эти результаты показывают эффективность предложенного подхода. На рис. 2 представлены рассчитанные волновые поля в разные моменты времени.
Рис. 2. Волновое поле в фиксированные моменты времени Благодарности
Работа была выполнена совместно с Московским Научноисследовательским центром компании Schlumberger и частично при поддержке грантов РФФИ 08-05-0265, 10-01-92604, 10-05-00233, 10-05-00337.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Berenger, J.P. A perfect matched layer for the absorption of electromagnetic waves // Journal of Comp. Physics. - 1994. - 114. - P. 185-200.
2. Chan, T. and Mathew, T.P. Domain decomposition // Acta Numerica. - 1994. - P. 61143.
3. Collino F., Ghanemi S. and Joly P. Domain Decomposition Method for Harmonic Wave Propagation // INRIA, Rapport de recherche. - 1998. - № 3473.
2. Collino F., Fouquet T. and Joly P. A Conservative Space-time Mesh Refinement Method for the 1-DWave Equation. Part 1: Construction // Numer. Math. - 2003. - vol. 95. - P. 197-221.
3. Joly P., Collino F. and Fouquet T. Analyse numerique d’une methode de raffinement de maillage espace-temps pour l’equation des ondes // INRIA, Rapport de recherche. - 1998. -3474.
4. Konyukh G.V. and Mikhailenko B.G. Application of the Laguerre Integral Transforms for solving Dynamic Seismic Problems // Bull. of the Novosibirsk Computing Center, Series: Mathematical Modeling in Geophys. - 1998. - № 5. - P. 71-92.
5. Martin J. Gander, Laurence Halpern and Frederic Nataf. Optimized Schwarz Methods // 12th International Conference on Domain Decomposition Methods. - 2001.
6. Plessix. A Helmholtz iterative solver for 3D seismic-imaging problems // Geophysics. - 2007.
7. Suetin P.K. Classical orthogonal polynomials // Moscow:Nauka. - 1974.
8. Гилбо Ж., Ланда Е., Решетова Г.В., Хайдуков В.Г., Чеверда В.А. Численное моделирование сейсмических волновых полей в двумерно-неоднородных упругих
разномасштабных средах (карстовые включения)// Технологии сейсморазведки. - 2008. -№3. - C. 19 - 28.
9. Virieux, J., Operto, S., Ben-Hadj-Ali, H., Brossier, R., Etienne, V., Sourbier, F., Giraud, L. and Haidar, A. Seismic wave modeling for seismic imaging // The Leading Edge 28(5). - 2009. - 538-544.
10. Решетова Г.В., Чеверда В.А. Использование преобразования Лагерра для построения идеально подходящих поглощающих слоев без расщепления // Математическое моделирование. - 2006. - Т. 18. - №10. - C. 91 - 101.
© М.А. Белоносов, 2010
