Разностные методы высокого порядка точности для решения акустического волнового уравнения и уравнений анизотропной упругости

Довгилович, Леонид Евгеньевич

Разностные методы высокого порядка точности для решения акустического волнового уравнения и уравнений анизотропной упругости тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Довгилович, Леонид Евгеньевич АВТОР

кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ

Москва МЕСТО ЗАЩИТЫ

2013 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

Диссертация по математике на тему «Разностные методы высокого порядка точности для решения акустического волнового уравнения и уравнений анизотропной упругости»

Автореферат диссертации на тему "Разностные методы высокого порядка точности для решения акустического волнового уравнения и уравнений анизотропной упругости"

005533748 „

На правах рукописи

Довгилович Леонид Евгеньевич

Разностные методы высокого порядка точности для решения акустического волнового уравнения и уравнений анизотропной упругости

01.01.07 - Вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

2 6 СЕН 2013

Москва-2013

005533748

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

Научный руководитель: Доктор физико-математических наук

Софронов Иван Львович

Официальные оппоненты: доктор физико-математических наук

Нечепуренко Юрий Михайлович Институт вычислительной математики Российской академии наук, ведущий научный сотрудник

кандидат физико-математических наук Хохлов Николай Игоревич кафедра информатики Московского фи-зикотехнического института, заведующий лабораторией

Ведущая организация: Институт прикладной математики им. М. В.

Келдыша РАН

Защита состоится « / 9"» Ъу^ -Я_ 2013 г. в _ часов

на заседании диссертационного совета Д 212.156.05 при Московском физико-техническом институте (государственном университете) по адресу: 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9, ауд. 903 КПМ.

С диссертацией можно ознакомиться в библиотеке Московского физико-технического института (государственного университета).

Автореферат разослан « »_сек %_2013 г.

Ученый секретарь диссертационного совета,

Федько О. С.

Общая характеристика работы

Актуальность темы исследования

В настоящее время в подавляющем большинстве случаев задачи моделирования сейсмических полей описываются скалярным волновым уравнением или системой линейных уравнений упругости. Специфика состоит в огромном объеме обрабатываемой и получаемой информации. Например, если говорить о сеточных методах, то размеры сеток могут достигать нескольких тысяч узлов в каждом направлении, а количество правых частей, т.е. отдельных задач, — десятков тысяч. Неудивительно, что такие расчеты требуют практически предельно возможных объемов оперативной памяти, и могут вестись неделями и месяцами на высокопроизводительных вычислительных системах. Типичными классами задач сейсмического моделирования с тысячами положений источников являются а) расчеты эталонных волновых полей и сейсмограмм — так называемая синтетика (synthetics); б) обратные задач сейсморазведки, где на каждой итерации решаются прямые задачи; в) задачи миграции в обратном времени. В последнее десятилетие резко возросла сложность моделируемых объектов как за счет учитываемых физических свойств (анизотропия, приповерхностные аномалии упругих свойств породы и т.п.), так и за счет геометрических особенностей (рельеф, сложная структура пластов и т.п.). Поэтому создание методов и алгоритмов, позволяющих решать такие задачи, является актуальной проблемой для вычислительной математики. Также актуальной является и возможность многократного уменьшения требуемой памяти и времени счета при использовании этих новых подходов там, где применяются традиционные.

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

Известны три основных класса методов, используемых при моделировании волновых процессов в неоднородных средах: метод конечных разностей (МКР), метод спектральных элементов (МСЭ), и спектральные методы (СМ). Указанные выше особенности 1), 2), 3) вполне определенно указывают на необходимость использования алгоритмов высокого порядка точности на криволинейных сетках для гиперболических задач с гладкими переменными коэффициентами. В применении к задачам сейсмики наиболее развитым в последнее десятилетие оказался МСЭ. Сочетание гибкости конечно-элементного представления сложной геометрии и спектральной аппроксимации уравнений внутри отдельного элемента привело к созданию достаточно универсальных и эффективных алгоритмов высокого порядка точности (эффективность алгоритма здесь означает использование как можно меньшего объема вычислительных ресурсов для достижения заданной точности решения). Построение аналогичных по гибкости алгоритмов на основе МКР высокого порядка точности для уравнений упругости является предметом данной работы. О СМ заметим лишь, что их применение в геофизике пока мало заметно.

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

Акустическое волновое уравнение рассматривается для частного случая ап-

проксимации оператора Лапласа на прямоугольных сетках с относительно простыми граничными условиями. Чтобы получить более эффективный алгоритм по сравнению с общепринятым центрально-разностным подходом, исследуются неявные центрально-разностные операторы (компактные схемы) вплоть до 20-го порядка точности для вычисления вторых производных.

Для системы Навье уравнений анизотропной упругости предлагается новый высокоточный конечно-разностный метод решения на криволинейных сетках. Его создание потребовало проведения предварительного анализа для выбора типа сетки — традиционно используемая разнесенная сетка или обычная. Оказалось, что обычная сетка позволяет экономить память при высоких порядках аппроксимации. Описываемый алгоритм обладает точностью вплоть до 8-го порядка во внутренних точках криволинейной сетки и четвертым порядком точности на границах. Интегрирование по времени реализовано явными схемами второго порядка точности, устойчивыми для любого невырожденного преобразования координат при соблюдении условия CFL (Courant-Friedrichs-Lewy). Построение схем основано на технологии суммирования по частям (Summation by Parts - SBP) и слабой формы учета граничных условий (Simultaneous Approximation Term - SAT). Для того, чтобы обеспечить максимальную скорость вычислений на современных высокопроизводительных системах, включая графические ускорители, алгоритм использует последовательное применение операторов первых производных на обычных (не разнесенных) сетках при аппроксимации операторов вторых и смешанных производных системы Навье уравнений анизотропной упругости. Избежать паразитных решений, свойственных центрально-разностным операторам на обычных сетках, позволяет аппроксимация первых производных разностями, сдвинутыми вперед и назад. Это также благоприятно сказалось на возможности аппроксимировать точечные источники, см. особенность 4) из вышеприведенного списка. Что касается особенности 5), то она учтена в рамках многоблочного подхода, см. [2], но этот вопрос не рассматривается в диссертации.

Отметим, что теория, развитая в работе при построении разностных схем для

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

Целью диссертационной работы является разработка эффективных конечно-разностных методов высокого порядка точности для решения акустического волнового уравнения и системы Навье нестационарных уравнений анизотропной упругости в трехмерных неоднородных средах.

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

Теоретическая и практическая значимость.

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

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

Реализованный параллельный алгоритм расчета сейсмических волн является основой комплекса программ миграции в обратном времени (Reverse Time

Migration - RTM) и обратной задачи сейсмики (Full Waveform Inversion - FWI), разработанного совместно с В.Г. Байдиным и И.Л. Софроновым. Комплекс применяется в Московском научно-исследовательском центре "Шлюмберже".

Положения, выносимые на защиту:

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

2. Построен параллельный алгоритм высокого порядка точности на криволинейных сетках для моделирования распространения волн в анизотропных упругих неоднородных средах. На основе многочисленных тестовых расчетов исследованы его точностные характеристики.

3. Предложен алгоритм разбиения на подобласти для параллельных вычислений с использованием неявных центрально-разностных операторов (компактных схем), основанный на перекрытии сеток. Теоретически и численно показано, что ширина перекрытия сопоставима с шириной шаблона.

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

Апробация результатов

Результаты работы были доложены, обсуждены и получили одобрение специалистов на следующих научных конференциях:

1. 52-ая, 53-ая и 55-ая конференции МФТИ, Москва-Долгопрудный, 2009, 2010, 2012;

2. XVIII и XIX Всероссийские конференции, посвященные памяти К.И. Бабен-ко, Абрау-Дюрсо, 2010, 2012;

3. 74-ая ежегодная конференция и выставка European Association Geoscientists & Engineers, Копенгаген, 2012;

4. Международная конференция «Разностные схемы и их приложения» посвященная 90-летию профессора B.C. Рябенького, Москва, 2013;

5. 75-ая ежегодная конференция и выставка European Association Geoscientists & Engineers, Лондон, 2013.

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

1. Институт вычислительной математики Российской академии наук, Москва, 2013;

2. Институт физики Земли им. О.Ю. Шмидта РАН, Москва, 2013;

3. Институт прикладной математики им. М.В. Келдыша РАН, Москва, 2013;

4. Московский физико-технический институт, кафедра информатики, Москва, 2013;

5. Московский научно-исследовательский центр "Шлюмберже", Москва, 2010-2013.

Публикации

Материалы диссертации опубликованы в 11 печатных работах, из них две [7, 8] - статьи в изданиях из списка, рекомендованного ВАК РФ. Оформлен патент [11].

Личный вклад автора.

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

Структура и объем диссертации.

Диссертация состоит из введения, шести глав, заключения, библиографии и приложения. Общий объем диссертации - 120 страниц. Библиография включает 96 наименований.

Содержание работы

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

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

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

Рассматривается одномерное волновое уравнение:

1 д2и д2и

где и(1, х) - неизвестная функция решения, А{х) > 0 - гладкая функция скорости распространения возмущений, х) - функция правой части. Для аппроксимации используется равномерная сетка по времени и на отрезке. Схема по времени — стандартная центрально-разностная второго порядка. А для приближения второй производной по пространству рассматриваются два семейства конеч-

но разностных операторов: явные центрально-разностные операторы и неявные центрально-разностные операторы

В параграфе 2.1 вводятся семейства рассматриваемых операторов дискретизации второй производной и ставится следующая задача: построение эффективного конечно-разностного алгоритма решения-уравнения (1) па равномерной сетке. Для этого формулируются критерии сравнения эффективности разных способов аппроксимации: 1/Точность решения при заданном объеме используемой памяти; 2)Объем используемой памяти, обеспечивающей заданную точность; 3)Ко-личество операций с плавающей точкой, затраченных на решение задачи с заданной точностью. На практике приходится находить компромисс для одновременного выполнения этих критериев, а также для учета эффективности программной реализации операторов, архитектуры вычислительной системы и т.п. Например, для минимального количества операций с плавающей точкой при малых точностях лучше использовать операторы с малым порядком аппроксимации. Но это влечет за собой использование очень мелкой сетки, что может привести к слишком большим объемам требуемой памяти для вычислений.

Рассматриваются два вида операторов аппроксимации второй производной. Первый вид - это центрально-разностные операторы (СБОр):

з=1

Второй вид - это неявные центрально-разностные операторы (1СООр), часто называемые компактными схемами:

(2)

В параграфе 2.2 подробно рассматривается периодический случай: проводятся асимптотический и спектральный анализ ошибок аппроксимаций. Показано, что использование 1С DO с трехдиагональными операторами в левой части предпочтительнее с точки зрения первых двух критериев. Из спектрального анализа заключаем, что оператор 1С DO 12 может считаться близким к оптимальному для практического применения при величинах погрешности вычисления второй производной е = 10~3 10~6.

Описание использования неявных центрально-разностных операторов для решения непериодических задач приводится в параграфе 2.3. В случае однородных условий Дирихле и Неймана предлагается воспользоваться симметрией и антисимметрией решения для вычислении значений оператора в приграничных точках соответственно. В общем случае граничных условий предлагается использовать значения вычисленные с помощью CDOp как условия Дирихле для граничных узлов отрезка. Доказывается теорема 2 о том, что ошибка, вносимая неправильной постановкой условий Дирихле для (3), убывает как а\ где |а| < 1, a г -расстояние до границы в точках.

В параграфе 2.4 приводится способ разбиения на подобласти для параллельных вычислений с использованием неявных центрально-разностных операторов. Благодаря сильному диагональному преобладанию обращаемого оператора и результатам теоремы 2 оказалось возможным применять те же принципы разбиения вычислительной области при реализации параллельных вычислений на системах с распределенной памятью, что и при использовании центрально-разностных операторов. В случае 1С DO область перекрытия соседних сеток ненамного больше, чем для С DO того же порядка (а с учетом, того, что для обеспечения заданной точности порядок берется меньшим чем у С DO, область перекрытия не вырастает). На практике при высоких порядках аппроксимации ширину перекрытия можно брать равной ширине шаблона правой части 1С DO.

Оценка количества операций с плавающей точкой для достижения заданной точности и сравнение рассматриваемых операторов с этой точки зрения приводят-

ся в параграфе 2.5. Исследован вопрос эффективного вычисления 1СЮО с трех-диагональными операторами в левой части. Продемонстрированы преимущества использования 1СЮО по третьему критерию. Обсуждения результатов и выводы данной главы приводятся в параграфе 2.6.

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

Рассматривается волновое уравнение на отрезке:

где м(£, х) - неизвестная функция решения, а(х) > 0 и Ъ{х) > 0 гладкие функции, /(£, х) - функция правой части.

Перечислим вначале ряд требований к дискретизации задачи (4), чтобы объяснить, почему понадобилось построение новых разностных схем. Во-первых, мы хотим построить консервативную разностную схему и тем самым обеспечить устойчивость счета на большие времена моделирования. Во-вторых, мы хотим аппроксимировать старшие производные по пространству путем последовательного применения первых производных (это позволит построить экономную реализацию соответствующего многомерного алгоритма на параллельных вычислительных системах). В третьих, мы не хотим использовать разнесенные сетки, чтобы избежать проблем с интерполяцией коэффициентов на полуцелые узлы. В четвертых, мы хотим, чтобы схема выдавала неосциллирующее решение в случае нарушения гладкости коэффициентов и/или в случае точечной правой части.

Взяв за основу идеологию БВР+БАТ, мы провели необходимое обобщение этого подхода, чтобы построить желаемое семейство схем высокого порядка ап-

(4)

проксимации по пространству. Особенностью схемы является последовательное использование разностных операторов производной 1-го порядка , £>+) для аппроксимации старших производных в уравнениях 2-го порядка.

В параграфе 3.1 обсуждаются существующие конечно-разностные способы решения уравнения (4).

Одна из главных проблем последовательной аппроксимации на обычных сетках анализируется для периодического случая в параграфе 3.2. Она состоит в том, что аппроксимация второй производной в волновом уравнении при помощи последовательного применения центрально-разностных операторов первой производной на обычных (не разнесенных) сетках приводит к проблеме чет-нечет, которая появляется при недостаточной гладкости решения и приводит к наличию нефизичных волн высокой частоты. Показано, что для всех рассмотренных порядков центрально-разностный оператор первой производной имеет нулевое собственное число при максимальной различимой гармонике, это и объясняет появление паразитных высокочастотных волн в решении. Используемое нами решение этой проблемы формулируется в параграфе 3.3. Мы применяем несимметричные операторы первой производной и т.е. в (4) мы аппроксимируем Л ~ 0~В0+. Несимметричность операторов заданного порядка Р достигается сдвигом шаблонов на одну точку, т.е. — £ + 1, £ + 1 для Вь и

у—тг — 1, j — lj для D~. Известно, что все собственные числа несимметричного оператора первой производной для высоких частот хорошо отделены от нуля, что обеспечивает отсутствие проблемы чет-нечет.

В параграфе 3.4 формулируется предложенное нами обобщение SBP подхода аппроксимации первой производной на отрезке в непериодическом случае, для сдвинутых разностных операторов из параграфа 3.3. Для этого на отрезке [0,1] вводится равномерная сетка из N узлов с шагом сетки h = . Аппроксимация первых производных из (4) будет проводится с помощью двух сеточных операторов D~ и D', удовлетворяющих условию SBP в следующем, модифицированном

нами, виде:

(й,0+г>)н + (0~й,у)н = «л- Уй, V е И,м, (5)

где Н = Нт > 0 - некоторая симметричная диагональная положительно определенная матрица. Мы вычислили коэффициенты операторов и В ' для условия (5) со значениями параметра Р = 4,6,8. Во внутренних точках отрезка операторы

01 и совпадают с операторами на сдвинутых шаблонах Р-то порядка, которые были описаны в параграфе 3.3. В Р точках возле границы отрезка (включая

саму границу) операторы £>+ и имеют порядок аппроксимации не меньше чем

2 •

Описание аппроксимации граничных условий Дирихле, Неймана и условий поглощения энергии на концах отрезка приводится в параграфе 3.5.

На основе разработанной пространственной аппроксимации в параграфе 3.6 строится явная по времени консервативная разностная схема. Для различных граничных условий в теоремах 3-6 доказана устойчивость и сходимость предложенных алгоритмов при выполнении условия СБЬ.

Обсуждение результатов и выводы данной главы приводятся в параграфе 3.7. В четвертой главе предлагается конечно-разностный метод решения трехмерного акустического волнового уравнения на прямоугольной сетке в параллелепипеде - одной из классических задач сейсмического моделирования. Рассматривается волновое уравнение в параллелепипеде:

1 д2и ^ д2и

= 9 (6)

г—1 1

где - неизвестная функция решения, А(х) > 0 - гладкая функция скоро-

сти распространения возмущений, д(£, х) — функция правой части. Для аппроксимации по времени используется стандартная явная схема второго порядка с постоянным шагом. Для аппроксимации оператора Лапласа используются неявные центрально-разностные операторы высокого порядка точности вычисления второй производной в каждом направлении.

В параграфе 4.1 формулируется задача данной главы следующего вида: построение эффективного конечно-разностного алгоритма решения уравнения (6) на равномерной сетке.

В параграфе 4.2 приводится дискретизация уравнения (6). Построенная схема обладает пространственным порядком вплоть до Р = 20. По времени используются явные схемы второго порядка. Особенностью схемы является применение неявных центрально-разностных операторов второй производной (компактных схем), что позволило увеличить эффективность по сравнению с общепринятым подходом на явных центрально-разностных операторах.

В параграфе 4.3 обсуждаются вопросы реализации. Для эффективной работы представленного алгоритма на высокопроизводительных системах с распределенной памятью используется распараллеливание путем разбиение на кубические подобласти. Это возможно благодаря утверждению теоремы 2 из главы 2.

В параграфе 4.4 приводятся численные эксперименты подтверждающие теоретический анализ предшествующих параграфов. На примере задачи расчета типовых трехмерных волновых полей для сейсмической миграции специализированным параллельным кодом с высокой степенью оптимизации получено, что минимальные затраты времени решения с использованием центрально-разностных схем достигаются при 20-ом порядке аппроксимации. Использование неявных центрально-разностных операторов для аппроксимации второй производной дает дополнительный выигрыш времени счета: в два раза для схем с 1СВО\ч по сравнению с С-ООго-

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

В пятой главе предлагается конечно-разностный метод решения трехмер-

ной системы Навье уравнений анизотропной упругости в неоднородных средах на криволинейных сетках. Аппроксимация проводится для волнового уравнения Навье в параметрических координатах в дивергентном виде, что позволяет получать криволинейные сетки в физических координатах, адаптированные к геометрическим и скоростным неоднородностям среды. Для этого разработан алгоритм на основе последовательной аппроксимации вторых и смешанных производных несимметричными разностными операторами первой производной высокого порядка точности на обычной (не разнесенной) сетке, построенными в главе 3, что позволяет избежать паразитных «чет-нечет» решений. Выписаны также аппроксимации условий Дирихле, условий свободной поверхности и условий поглощения энергии для произвольной гладкой границы и анизотропной гладкой среды. Для аппроксимации по времени используется стандартная явная схема второго порядка с постоянным шагом. Доказана ее устойчивость и сходимость при соблюдении условия Куранта-Фридрихса-Леви.

В параграфе 5.1 перечисляются основные конечно-разностные подходы к задачам сейсмики. В параграфе 5.2 формулируется задача динамической упругости и исследуется вопрос выбора конечно-разностной аппроксимации.

В данном параграфе продемонстрированы основные преимущества разработанного подхода над существующими конечно-разностными алгоритмами, применяемыми для расчетов анизотропных сейсмических полей. Так алгоритмы на разнесенных сетках Лебедева требуют большего количества вычислительных ресурсов для достижения той же точности. Кроме того, вопрос об аппроксимации граничных условий с высоким порядком точности на разнесенных сетках остается открытым. Консервативный алгоритм на обычных сетках с использованием БВР подхода в изотропном случае уже описан в литературе. Он имеет порядок Р — 4 (имеется ввиду порядок во внутренних точках). Для изотропного случая и криволинейной границы данный метод был представлен только с порядком Р = 2. В рамках этого подхода невозможно построить «разумный» алгоритм порядка выше чем Р = 6, из-за слишком маленького шага по времени необходимого для Р = 8.

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

В параграфе 5.3 приводится схема пространственной дискретизации. Используется трехмерная равномерная сетка в параметрических координатах 0 < < 1 с Агг, i = 1, 2,3 точек в каждом направлении. Для аппроксимации первых производных используются пары трехмерных операторов и £)~, которые строятся как тензорное произведение одномерных БВР операторов из главы 3.

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

= - ^н-1 + ¿,¿ = 1,2,3,

<4 = Су,¿.¿л* = 1,2,3, (7)

е% = \ (ткзО+и'1 + ТпБ+и^ , ¿,.7 = 1,2,3,

где у!}, а^, - сеточные функции, аппроксимирующие щ, ггу, соответственно; I, Р — диагональные операторы со значениями якобиана преобразования и плотности соответственно, Ту — диагональные операторы со значениями Сг]ц -диагональные операторы со значениями коэффициентов закона Тука, Н = Н1 Я20//\ //' - матрицы скалярных произведения из (5), а б}/,- трехмерный оператор, который строится как тензорное произведение единичных одномерных операторов I и соответствующих одномерных операторов (} = (Иад(—1,0,0,..., 0, 0,1).

В параграфе 5.5 к (7) применена явная схема аппроксимации по времени второго порядка и в теоремах 7, 8 доказана устойчивость и сходимость полученных алгоритмов при выполнении условия СИЬ.

Вопросы связанные с параллельной реализацией алгоритма, зависимости размера необходимой памяти от типа симметрии в законе Гука, способы экономии памяти и т.п. затронуты в параграфе 5.6. Обсуждение результатов и выводы

данной главы приводятся в параграфе 5.1.

В шестой главе приведены примеры вычислительных экспериментов с использованием конечно-разностной схемы из главы 5.

В параграфе 6.1 рассматривается задача с гладкой анизотропной средой в области с криволинейными границами и условиями свободной поверхности на них, с целью показать сеточную сходимость и устойчивость. Численно подтверждено, что построенный алгоритм сходится с 4-м порядком точности во всей вычислительной области вплоть до границы. При рассмотрении сходимости в области, отделенной от границы, мы наблюдаем увеличение порядка сходимости вплоть до 8-го на некотором времени моделирования, что согласуется с фактом использования операторов 8-го порядка точности во внутренних точках.

В параграфе 6.2 приводятся результаты по моделированию двумерной задачи Лэмба, имеющей аналитическое решение. Продемонстрирована сходимость к этому решению, в том числе и при использовании криволинейных координат (см. рис. 1Ь). Показано преимущество использования операторов высокого порядка. Также были исследованы свойства разностной схемы при моделировании точечного по пространству источника одной точкой сетки и гладким сигналом по времени. Показано, что в отдалении от точки источника на расстояние порядка ширины шаблона решение сходится с заданным порядком (4-ым в данном случае), см. рис. 1а.

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

18 14

ИО 6

•(1=12 И=6

12 60 132 300 420 540 660

(а)

Рис. 1. Задача Лэмба: (а) - погрешность решения на поверхности в зависимости от расстояния от точечного источника, шаги сетки 12 м и 6 м; (Ь) — модуль вектора решения для случая наклонной поверхности в моменты времени 0.44 с;

(а)

(Ь)

Рис. 2. Геометрия модельной задачи: (а) - поверхность; (Ь) - двумерный срез расчетной сетки;

Обсуждение результатов и выводы данной главы мы приводим в параграфе 6.4.

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

(а)

Рис. 3. Модуль вектора решения: (а) - изотропный случай; (Ь) - анизотропный случай;

Основные результаты работы

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

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

3. Проведен анализ конечно-разностных аппроксимаций высокого порядка на разнесенных и обычных сетках с целью минимизации вычислительных ресурсов, необходимых для решения волнового уравнения Навье в криволинейных координатах для анизотропных сред. Показано, что использование обычных сеток предпочтительнее для схем высокого порядка (> 04).

4. Предложена конечно-разностная схема высокого порядка точности для решения трехмерных задач динамической упругости в анизотропных неоднородных средах на криволинейных сетках, адаптированных к геометри-

(Ь)

ческим и скоростным неоднородностям среды. Схема обладает пространственным порядком вплоть до 08 во внутренних точках и вплоть до 04 на границах. Проанализированы устойчивость и сходимость.

5. С целью минимизации числа операций на точку и эффективного распараллеливания на многоядерных вычислительных системах использован подход последовательного применения разностных операторов первого порядка для аппроксимации системы уравнений второго порядка. Для устранения паразитных «чет-нечет»-решений разработана схема аппроксимации на сдвинутых вперед и назад шаблонах.

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

7. Реализован эффективный параллельный алгоритм высокого порядка аппроксимации для решения трехмерного волнового уравнения Навье на криволинейных сетках в неоднородных анизотропных средах.

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

9. Реализованы алгоритмы решения прямых задач для комплекса программ миграции в обратном времени, разработанного совместно с В.Г. Байдиным и И.Л. Софроновым.

Публикации автора по теме диссертации

1. Dovgilovich L. High-order FD Method on Curvilinear Grids for Anisotropic Elastodynamic Simulations // 75th EAGE Conference & Exhibition. —

2013.— P. 1—5.— URL: http://www.earthdoc.org/publication/ publicationdetails/?publication=68514.

2. Sofronov IL, Zaitsev NA, Dovgilovich L et al. Multi-block FD Method for 3D Geophysical Simulation with Explicit Representation of Sub-horizontal Interfaces // 74th EAGE Conference & Exhibition.— 2012,— P. 1-5.— URL: http://www.earthdoc.org/publication/publicationdetails/ ?publication=58892.

3. Довгилович JI.E. Анализ явных и неявных конечно-разностных операторов в применении к задачам распространения волн в анизотропной упругой среде // Тезисы докладов XIX Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», посвященной памяти К.И. Бабенко.— Москва: Институт прикладной математики им. М.В. Келдыша, 2012. — С. 36-37.

4. Довгилович JI.E. Параллельный конечно-разностный алгоритм моделирования распространения упругих волн в трехмерной анизотропной гетерогенной среде с условиями свободной поверхности на криволинейных сетках // Современные проблемы фундаментальных и прикладных наук. Часть VII. Управление и прикладная математика: Труды 55-й научной конференции МФТИ.— Т. 2. — Москва: МФТИ, 2012, — С. 87-88.

5. Довгилович JI.E. Анализ производительности алгоритма решении обратной задачи сейсморазведки методом миграции в обратном времени с применением компактных схем // Современные проблемы фундаментальных и прикладных наук. Часть VII. Управление и прикладная математика: Труды 53-й научной конференции МФТИ. — Т. 3, — Москва: МФТИ, 2010. — С. 19-20.

6. Довгилович Л.Е. О решении обратной задачи сейсморазведки методом миграции в обратном времени с применением компактных схем // Современные проблемы фундаментальных и прикладных наук. Часть VII. Управление и

прикладная математика: Труды 52-й научной конференции МФТИ.— Т. 3.— Москва: МФТИ, 2009, — С. 122-124.

7. Довгилович JI.E, Софронов И.Л. Анализ явных и неявных центрально-разностных операторов для вычисления второй производной на равномерных сетках // Труды МФТИ, — 2013. — Т. 5, № 2, — С. 175-182.

8. Довгилович Л.Е, Софронов И.Л. Конечно-разностный метод высокого порядка точности расчета волновых полей в анизотропных средах // Технологии Сейсморазведки. — 2013.— № 2. — С. 24—30.

9. Довгилович Л.Е, Софронов И.Л. Аппроксимация вторых и смешанных производных SBP подходом на основе несимметричных первых разностей // Сборник научных трудов Международной конференции «Разностные схемы и их приложения» посвященной 90-летию профессора B.C. Рябенького. — Москва: Институт прикладной математики им. М.В. Келдыша, 2013.— С. 52-53.

10. Довгилович Л.Е, Софронов И.Л. О трехуровневом параллельном алгоритме обратной миграции во временной области (RTM) на основе компактных схем // Тезисы докладов XVIII Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», посвященной памяти К.И. Бабенко.— Москва: Институт прикладной математики им. М.В. Келдыша, 2010. — С. 31-32.

Патенты:

11. Baydin Vasily, Dovgilovich Leonid, Lin Kui et al. Reverse time migration model dip-guided imaging. Appl. No.: 13/725,154; Pub. No.: US 2013/0182538 Al; за-явл. 21.12.2012; опубл. 18.07.2013, 14c.

Довгилович Леонид Евгеньевич

Разностные методы высокого порядка точности для решения акустического волнового уравнения и уравнений анизотропной

упругости

Автореферат

Подписано в печать 12.09.2013. Формат 60 х 841/16. Усл. печ. л. 1,0. Тираж 100 экз. Заказ № 266.

Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный

университет)» Отдел оперативной полиграфии «Физтех-полиграф» 141700, Московская обл., г. Долгопрудный, Институтский пер., 9.

Текст научной работы диссертации и автореферата по математике, кандидата физико-математических наук, Довгилович, Леонид Евгеньевич, Москва

Московский физико-технический институт (государственный университет) Кафедра вычислительной математики

На правах рукописи УДК 519.63

04201361484

ДОВГИЛОВИЧЛеонид Евгеньевич

Специальность 01.01.07 - Вычислительная математика

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель

доктор физико-математических наук

Софронов Иван Львович

Москва-2013

Содержание

Введение..............................................................................................5

Глава 1. Моделирование волновых полей в задачах геофизики............................10

1.1. О роли моделирования в геофизике......................................................10

1.2. Рассматриваемые модели ................................................................12

1.3. Обзор методов..............................................................................13

Глава 2. Аппроксимация второй производной неявными центрально-разностными

операторами ....................................................................................16

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

2.1.1. Исследуемые операторы второй производной ................................17

2.2. Периодический случай....................................................................19

2.2.1. Асимптотический анализ........................................................19

2.2.2. Спектральный анализ............................................................22

2.3. Непериодический случай ................................................................26

2.3.1. Однородные условия Дирихле и Неймана......................................26

2.3.2. Общий случай....................................................................27

2.4. Использование неявных центрально-разностных операторов для параллельных вычислений................................................................................30

2.5. Оценки количества операций с плавающей точкой для получения заданной точности 32

2.6. Обсуждение результатов..................................................................36

Глава 3. Консервативные разностные схемы высокого порядка точности по пространству для одномерного волнового уравнения................................................38

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

3.2. Проблема чет-нечет........................................................................40

3.3. Сдвинутые шаблоны......................................................................44

3.4. 8ВР подход для операторов на сдвинутых шаблонах..................................46

3.5. Граничные условия........................................................................50

3.5.1. Условия Дирихле ................................................................50

3.5.2. Условия Неймана................................................................51

3.5.3. Условия Неймана с включением граничных точек............................52

3.5.4. Условия третьего рода и поглощения энергии ................................54

3.6. Консервативная разностная схема, сходимость и устойчивость......................56

3.7. Обсуждение результатов..................................................................59

Глава 4. Разностный метод высокого порядка для решения трехмерного волнового

уравнения на равномерной сетке ............................................................60

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

4.2. Дискретизация ............................................................................61

4.2.1. Граничные условия..............................................................62

4.3. Реализация..................................................................................63

4.4. Численные эксперименты................................................................63

4.4.1. Тестирование основных характеристик алгоритма. Отладка программы . . 63

4.4.2. Пример трехмерного расчета различными операторами......................64

4.5. Обсуждение результатов..................................................................65

Глава 5. Разностный метод высокого порядка точности для решения уравнений анизотропной упругости в криволинейных координатах........................................67

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

5.2. Выбор типа конечно-разностной аппроксимации для решения задачи................69

5.2.1. Консервативная форма записи уравнений анизотропной упругости в криволинейных координатах........................................................69

5.2.2. Последовательная и смешанная аппроксимация..............................70

5.2.3. Разнесенные и не разнесенные сетки..........................................71

5.3. Дискретизация уравнений анизотропной упругости в криволинейных координатах 75

5.4. Граничные условия........................................................................77

5.4.1. Условия Дирихле ................................................................77

5.4.2. Условия свободной поверхности................................................79

5.4.3. Условия поглощения энергии....................................................81

5.4.4. Комбинирование граничных условий..........................................83

5.5. Схема по времени, устойчивость и сходимость предложенного метода..............84

5.6. Реализация..................................................................................85

5.7. Обсуждение результатов..................................................................86

Глава 6. Численные эксперименты на задачах упругости..................................88

6.1. Сеточная сходимость ....................................................................88

6.2. Задача Лэмба..............................................................................94

6.3. Сейсмические модельные эксперименты................................................97

6.4. Обсуждение результатов.................................107

Заключение.............................................108

Список использованных источников..............................110

Приложение А. Значения коэффициент для операторов второй производной .....119

Введение

Актуальность темы исследования

В настоящее время в подавляющем большинстве случаев задачи моделирования сейсмических полей описываются скалярным волновым уравнением или системой линейных уравнений упругости. Специфика состоит в огромном объеме обрабатываемой и получаемой информации. Например, если говорить о сеточных методах, то размеры сеток могут достигать нескольких тысяч узлов в каждом направлении, а количество правых частей, т.е. отдельных задач, - десятков тысяч. Неудивительно, что такие расчеты требуют практически предельно возможных объемов оперативной памяти, и могут вестись неделями и месяцами на высокопроизводительных вычислительных системах. Типичными классами задач сейсмического моделирования с тысячами положений источников являются а) расчеты эталонных волновых полей и сейсмограмм - так называемая синтетика (synthetics); б) обратные задач сейсморазведки, где на каждой итерации решаются прямые задачи; в) задачи миграции в обратном времени. В последнее десятилетие резко возросла сложность моделируемых объектов как за счет учитываемых физических свойств (анизотропия, приповерхностные аномалии упругих свойств породы и т.п.), так и за счет геометрических особенностей (рельеф, сложная структура пластов и т.п.). Поэтому создание методов и алгоритмов, позволяющих решать такие задачи, является актуальной проблемой для вычислительной математики. Также актуальной является и возможность многократного уменьшения требуемой памяти и времени счета при использовании этих новых подходов там, где применяются традиционные.

Основные вычислительные особенности, связанные с рассматриваемыми классами задач, состоят в следующем: 1) расчет волновых полей должен вестись с контролируемой и достаточно высокой точностью; 2) допустимо использование достаточно гладких параметров среды при решении обратных задач и задач миграции в обратном времени (сильные градиенты в местах предполагаемых границ пластов при этом вполне возможны); 3) должен учитываться рельеф земной поверхности; 4) источники могут располагаться вблизи целевой зоны сейсморазведки, например, для систем наблюдений в скважинной сейсмике (речь идет о сопряженных уравнениях в теории обратных задач); 5) должна иметься возможность учета разрывов свойств среды (разломы, солевые тела, пласты и т.п.).

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

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

Акустическое волновое уравнение рассматривается для частного случая аппроксимации оператора Лапласа на прямоугольных сетках с относительно простыми граничными условиями. Чтобы получить более эффективный алгоритм по сравнению с общепринятым центрально-разностным подходом, исследуются неявные центрально-разностные операторы (компактные схемы) вплоть до 20-го порядка точности для вычисления вторых производных.

Для системы Навье уравнений анизотропной упругости предлагается новый высокоточный конечно-разностный метод решения на криволинейных сетках. Его создание потребовало проведения предварительного анализа для выбора типа сетки - традиционно используемая разнесенная сетка или обычная. Оказалось, что обычная сетка позволяет экономить память при высоких порядках аппроксимации. Описываемый алгоритм обладает точностью вплоть до 8-го порядка во внутренних точках криволинейной сетки и четвертым порядком точности на границах. Интегрирование по времени реализовано явными схемами второго порядка точности, устойчивыми для любого невырожденного преобразования координат при соблюдении условия CFL (Courant-Friedrichs-Lewy). Построение схем основано на технологии суммирования по частям (Summation by Parts - SBP) и слабой формы учета граничных условий (Simultaneous Approximation Term -SAT). Для того, чтобы обеспечить максимальную скорость вычислений на современных высокопроизводительных системах, включая графические ускорители, алгоритм использует последовательное применение операторов первых производных на обычных (не разнесенных) сетках при аппроксимации операторов вторых и смешанных производных системы Навье уравнений анизотропной упругости. Избежать паразитных решений, свойственных центрально-разностным one-

раторам на обычных сетках, позволяет аппроксимация первых производных разностями, сдвинутыми вперед и назад. Это также благоприятно сказалось на возможности аппроксимировать точечные источники, см. особенность 4) из вышеприведенного списка. Что касается особенности 5), то она учтена в рамках многоблочного подхода, см. [53], но этот вопрос не рассматривается в диссертации.

Отметим, что теория, развитая в работе при построении разностных схем для системы Навье уравнений анизотропной упругости, может быть использована и в модели акустического волнового уравнения, если, например, возникнет необходимость применения криволинейных сеток для учета топографии поверхности Земли.

Теоретическая и практическая значимость.

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

Реализованный параллельный алгоритм расчета сейсмических волн является основой комплекса программ миграции в обратном времени (Reverse Time Migration - RTM) и обратной задачи сейсмики (Full Waveform Inversion - FWI), разработанного совместно с В.Г. Байдиным и И.Л. Софроновым. Комплекс применяется в Московском научно-исследовательском центре "Шлюм-берже".

Положения, выносимые на защиту:

Апробация результатов

1. 52-ая, 53-ая и 55-ая конференции МФТИ, Москва-Долгопрудный, 2009, 2010, 2012;

2. XVIII и XIX Всероссийские конференции, посвященные памяти К.И. Бабенко, Абрау-Дюрсо, 2010, 2012;

3. 74-ая ежегодная конференция и выставка European Association Geoscientists & Engineers, Копенгаген, 2012;

5. 75-ая ежегодная конференция и выставка European Association Geoscientists & Engineers, Лондон, 2013.

1. Институт вычислительной математики Российской академии наук, Москва, 2013;

2. Институт физики Земли им. О.Ю. Шмидта РАН, Москва, 2013;

3. Институт прикладной математики им. М.В. Келдыша РАН, Москва, 2013;

4. Московский физико-технический институт, кафедра информатики, Москва, 2013;

5. Московский научно-исследовательский центр "Шлюмберже", Москва, 2010-2013. Публикации

Материалы диссертации опубликованы в 11 печатных работах, из них две [70, 71 ] - статьи в изданиях из списка, рекомендованного ВАК РФ. Подана и опубликована заявка на патент [96].

Личный вклад автора.

Структура и объем диссертации.

Глава 1. Моделирование волновых полей в задачах геофизики 1.1. О роли моделирования в геофизике

Сейсморазведка является наиболее распространенным поверхностным метод�