CC BY
44
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВОЛНОВЫХ ПРОЦЕССОВ В ТРЁХМЕРНЫХ УПРУГИХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ ПРЕОБРАЗОВАНИЯ ЛАГЕРРА
Михаил Андреевич Белоносов
Федеральное государственное бюджетное учреждение науки Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, проспект Академика Коптюга 3, инженер, e-mail: belonosovma@ipgg.sbras.ru
Сергей Александрович Соловьев
Федеральное государственное бюджетное учреждение науки Институт нефтегазовой геологии и геофизики Сибирского отделения РАН, 630090, Россия, г. Новосибирск, проспект Академика Коптюга 3, к.ф.-м.н., ведущий инженер, e-mail: solovevsa@ipgg.sbras.ru
В настоящей работе для численного моделирования сейсмических волновых полей в упругих неоднородных средах используется подход, основанный на применении интегрального преобразования Лагерра с последующей декомпозицией области. Благодаря тому, что преобразование Лагерра приводит к системе с отрицательно определённым эллиптическим оператором, независящим от параметра разделения, параллельные вычисления можно организовать с помощью аддитивного метода Шварца, а решение систем линейных алгебраических уравнений в каждой подобласти с помощью LU-разложения. В работе исследуются возможности применения данного подхода в трёхмерном случае.
Ключевые слова: преобразование Лагерра, аддитивный метод Шварца, волновые поля, параллельные вычисления.
NUMERICAL SIMULATION OF WAVE PROPAGATION IN 3D ELASTIC MEDIUM WITH LAGUERRE TRANSFORM
Mikhail A. Belonosov
A.A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3, Akademika Koptyuga Prosp., engineer, e-mail: belonosovma@ipgg.sbras.ru
Sergey A. Solovev
A.A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3, Akademika Koptyuga Prosp., PhD in physics and mathematics, senior engineer, email: solovevsa@ipgg.sbras.ru
In this paper, we apply an approach for numerical simulation of elastic waves in heterogeneous medium, based on the implementation of integral Laguerre transform with subsequent domain decomposition. Following the Laguerre transform, we obtain a system with strictly negative definite elliptic operator, which doesn’t depend on separation parameter. Therefore parallel calculat ions can be organized by means of the additive Schwarz method and systems of linear algebraic equations in each subdomain can be solved by means of LU factorization. In the paper we study this approach in 3D case.
Key words: Laguerre transform, additive Schwarz method, wavefield, parallel calculations.
Введение
Проведение полномасштабного численного моделирования сейсмических волн в реалистичных трёхмернонеоднородных средах невозможно без использования кластерных систем с распределённой памятью. В настоящее время наиболее распространённым направлением для этого является использование явных конечно-разностных схем. Однако этот подход имеет ряд существенных недостатков, среди которых наиболее существенным является обмен данными между соседними процессорными элементами на каждом шаге по времени. Поэтому значительные усилия применяются для развития методов численного моделирования, основанных на отделении времени. Основное внимание здесь сосредоточено на отделении времени путём применения преобразования Фурье [1]. Однако при этом возникают существенные трудности при решении задач в реалистичных постановках и связаны они в основном со знаконеопределённо-стью пространственного оператора. Именно из-за этого скорость итерационных процессов решения систем линейных алгебраических уравнений (СЛАУ) неприемлемо низка.
В настоящей работе отделение времени проделывается путем применения интегрального преобразования Лагерра. Оно приводит к знакоопределённому эллиптическому оператору [2], что создаёт основу для возможного применения метода декомпозиции области и альтернирования по Шварцу [3] для решения полученной СЛАУ [4]. На этой базе разработан параллельный алгоритм численного моделирования сейсмических волновых полей и их свойства исследованы в двумерном случае [5,6]. Здесь же этот алгоритм применяется для трёхмерной системы уравнений динамической теории упругости.
Постановка задачи и алгоритм её решения
В некотором прямоугольном параллелепипеде рассмотрим систему уравнений динамической теории упругости для источника типа центр объёмного расширения с нулевыми начальными данными:
где р - плотность, С i jki - тензор жёсткости, х0 - координаты источника, f ( t) -излучаемый им импульс. Для того, чтобы избежать отражений от границы расчётной области она окаймляется идеально подходящим поглощающим слоем [7] или PML (аббревиатура от английского Perfectly Matched Layer).
Применение интегрального преобразования Лагерра по времени к системе (1) трансформирует её в систему
L ип = д п,
с оператором , обладающим следующими свойствами:
- строго отрицательно определённый;
- не зависит от параметра разделения .
Для дальнейшего решения задачи расчётная область разбивается на несколько перекрывающихся подобластей Рис. 1а и применяется аддитивный метод Шварца. Его сходимость гарантируется знакоопределённостью оператора. Основная идея этого метода состоит в том, что решение вычисляется не во всей расчётной области, а в каждой отдельной элементарной подобласти её пространственной декомпозиции с последующей организацией обменов между соседями. В частности, для отыскания решения в области , вводится разбиение последней на перекрывающиеся подобласти Ог и О2 (см. Рис. 1б). На первой итерации метода Шварца решение ищется в этих подобластях с произвольными граничными условиями Дирихле на границах и . На каждой последующей итерации на каждой из этих границ ставится условие, равное значению уже найденного решения задачи на предыдущей итерации из соседней подобласти.
Рис. 1. а. Декомпозиция трёхмерной расчётной области; б. Декомпозиция расчетной области с перекрытием для организации
альтернирований по Шварцу
Организация параллельных вычислений
Второе свойство системы (2) открывает возможность использования прямых методов решения СЛАУ, возникающих после её конечно -разностной аппроксимации. Таким образом мы имеем два вычислительных процесса:
1) решение СЛАУ прямыми методами на многоядерных кластерных узлах, используя OpenMP;
2) организация обменов между подобластями при выполнении альтернирований по Шварцу с использованием библиотеки MPI.
Для решения СЛАУ на каждом кластерном узле используется компонента PARDISO из математической библиотеки Intel®MKL. Процесс факторизации в данной компоненте распараллелен и оптимизирован для работы с многоядерными вычислительными системами с общей памятью (OpenMP).
В результате использования формата CSR (compress sparse row), предназначенного для эффективного хранения разреженных матриц, оказалось, что основной объём оперативной памяти требуется для хранения LU-разложения матрицы. Так в двумерном случае для решения задачи в области 1 6 О О X 1 6 О О точек нужно около 8 ГБ оперативной памяти. В трёхмерном же случае, для областей совсем небольших размеров требуется огромный объём оперативной
памяти. Так для 8 Гб мы можем использовать подобласти размера лишь 70x70x70 точек.
Процесс LU-разложения матриц, получаемых в результате аппроксимации трехмерных задач представляет собой отдельный вопрос, т.к. заполняемость LU-факторов ненулевыми элементами не может быть значительно уменьшена, даже с помощью известных пакетов как METIS либо SCOTCH. Данные пакеты используют наиболее эффективный алгоритм переупорядочивания строки столбцов исходной матрицы Nested Dissection (ND) [8], однако для трехмерных задач необходимо применять дополнительные подходы. Один из них основан на использовании свойств малого ранга, которыми обладают внедиагональные блоки LU-факторов, а также внедиагональные блоки дополнения Шура [9], что позволяет существенно сократить затраты на вычислительные ресурсы [10]. Оценка числа арифметических операций для факторизации системы, а также оценка затрат по памяти представлена в табл. 1.
Таблица 1
Затраты на вычисление LU-разложения матрицы при решении трехмерной задачи в области точек
LU-разложение с ND- алгоритмом с использованием свойств малого ранга
Число арифметических операций 0(N8) 0(N6) 0(N4)
Необходимый объем оперативной памяти 0(N5) 0(N4) 0(N3)
Численные эксперименты
Для того, чтобы оценить скорость сходимости аддитивного метода Шварца, была проделана серия численных экспериментов. Рассматривалась однородная среда, которая разбивалась на различное число подобластей как показано на Рис. 2. Для каждой такой декомпозиции решалась система (1) и замерялось число итераций в методе Шварца. В Табл. 2 представлен результат этих расчётов.
Рис. 2. Одномерная декомпозиция расчётной области
Таблица 2
Скорость сходимости метода Шварца
Число подобластей 2 3 4 5 6 7 8
Число итераций 3 4 5 6 7 8 9
Для слоистой среды 8 0 0 x 8 0 0 x 8 0 0 м, представленной на Рис. 3а (в
у„
верхнем слое V, = 3 60 0 м/с, в нижнем слое V, = 5 0 0 0 м/с; = ^=, р = 1 5 0 0
кг/;и3) были рассчитаны волновые поля. Область покрывалась равномерной сеткой с шагом 5 м. Источник типа центр объемного расширения, излучающий импульс Рикера с доминирующей частотой 30 Гц, расположен в центре этой области. Она разбивается на 27 перекрывающихся подобластей (3x3x3 ) (см. Рис. 1а). Каждая подобласть загружается на свой вычислительный узел, имеющий 8 Гб оперативной памяти. Результат моделирования в виде моментального снимка волнового поля в двух проекциях приведен на Рис. 3б.
Рис. 3. а. Вертикально-неоднородная слоистая среда; б. Результаты численного моделирования: моментальный снимок волнового поля
Численные эксперименты выполнялись на вычислительных системах Московского государственного университета и Межведомственного суперкомпью-терного центра.
Благодарности
Данная работа выполнена при поддержке грантов РФФИ 11 -05-00947, 1205-31008, 13-05-00076 и интеграционных проектов СО РАН № 127 и 130.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Plessix R.-E. A Helmholtz iterative solver for 3D seismic-imaging problems // Geophysics, 72. - 2007. - P. 185-194.
2. Mikhailenko B.G., Mikhailov A.A., Reshetova G.V. Numerical viscoelastic modeling by the spectral Laguerre method // Geophysical Prospecting. - 2003. - V. 51. - P. 37-48.
3. Мацокин А.М., Непомнящих С.В. Метод альтернирования Шварца в подпространствах // Изв. высших учебных заведений. - 1985. - Т. 29. - № 10. - С. 61-66.
4. Chan T., Mathew T.P. Domain decomposition // Acta Numerica. - 1994. - V. 3. - P. 61143.
5. Белоносов М.А., Костов К., Решетова Г.В., Соловьёв С.А., Чеверда В.А. Организация параллельных вычислений для моделирования сейсмических волн с использованием аддитивного метода Шварца // Вычислительные методы и программирование. - 2012. - Т. 13. -С. 525-535.
6. Belonosov M.A., Kostov C., Reshetova G.V., Soloviev S.A., Tcheverda V.A. Parallel numerical simulation of seismic waves propagation with Intel Math Kernel Library // Lecture Notes in Computer Science. Applied Parallel and Scientific Computing. - 2013. - V. 7782. - P. 153-167.
7. Решетова Г.В., Чеверда В.А. Использование преобразования Лагерра для построения идеально подходящих поглощающих слоев без расщепления // Математическое моделирование. - 2006. - Т. 18. - №10. - C. 91 - 101.
8. J. A. George, Nested dissection of a regular finite element mesh, SIAM J. Numer. Anal., 10 (1973), pp. 345-363.
9. Chandrasekaran, S., Dewilde, P., Gu, M., Somasunderam, N.: On the numerical rank of the off-diagonal blocks of Schur complements of discretized elliptic PDEs. SIAM J. Matrix Anal. Appl. 31, 2261-2290 (2010)
10. J. Xia, Robust and efficient multifrontal solver for large discretized PDEs, High-Perform. Sci. Comput., M. W. Berry et al. (eds.), Springer (2012), pp. 199-217.
© М.А. Белоносов, С.А. Соловьев, 2013
