Содержание к диссертации
Введение
Глава 1 Постановка задачи и ее аппроксимация 26
1.1 Постановка задачи 27
1.2 Дискретизация по пространству 31
1.3 Полудискретная система 35
1.4 Дискретизация по времени. Схема с весами 38
Глава 2 Схемы расщепления для вектора теплового потока 44
2.1 Схемы расщепления с треугольной факторизацией 46
2.1.1 Схема на основе попеременно-треугольной факторизации 47
2.1.2 Схема на основе факторизации типа SSOR 2.2 Схемы расщепления на основе алгоритма Удзавы 51
2.3 Схемы расщепления на основе скалярных схем расщепления для дивергенции теплового потока
2.3.1 Двумерный случай 55
2.3.2 Трехмерный случай 60
2.3.3 Аппроксимация потоковых схем расщепления на основе схем для дивергенции теплового потока 64
2.4 Устойчивость схем расщепления. Априорные оценки 69
2.4.1 Устойчивость в полунорме для потоковых схем расщепления. Общие результаты 70
2.4.2 Априорные оценки глобальной устойчивости для потоковой схемы расщепления на основе схемы переменных направлений 72
2.4.3 Априорные оценки глобальной устойчивости для потоковой схемы расщепления на основе схемы Дугласа-Гана 78
Глава 3 Вычислительные эксперименты 81
3.1 Сравнение схем по точности 82
3.1.1 Двумерный случай 83
3.1.2 Трехмерный случай
3.2 Априорные оценки и гладкость начальных данных 99
3.3 Сравнительный анализ производительности параллельных реализаций с MPI и ОрепМР 1 3.3.1 Описание MPI реализации 106
3.3.2 Описание MPI/OpenMP реализации с "почтальонами" 107
3.3.3 Сравнительный анализ параллельных реализаций: MPI, "простой" MPI/OpenMP реализации и MPI/ ОрепМР с "почтальонами" 111
Глава 4 Численное моделирование термохронологии аккре ционно-коллизионных процессов в литосфере Земли 117
4.1 Термохронологические реконструкции для неопротерозойской активной континентальной окраины
Сибирского кратона 119
4.1.1 Физическая модель 119
4.1.2 Математическая модель 124
4.1.3 Результаты моделирования 127
4.2 Моделирование термохронологии косой позднепалеозойской коллизии в зоне Таймыр Северноземельской складчатости 129
4.2.1 Физическая модель 130
4.2.2 Математическая модель 135
4.2.3 Результаты моделирования 143
Заключение 148
Литература
- Полудискретная система
- Схемы расщепления на основе скалярных схем расщепления для дивергенции теплового потока
- Априорные оценки и гладкость начальных данных
- Моделирование термохронологии косой позднепалеозойской коллизии в зоне Таймыр Северноземельской складчатости
Полудискретная система
Очевидно, билинейные формы a(-, ) и m(-, ) являются симметричными, положительно определенными на пространствах V и Q соответственно. Тогда обобщенная (по пространственным переменным) постановка задачи может быть сформулирована следующим образом: Найти (w,T) Є C(0,tf;V) х Cl(Q,tf;Q) такие, что
Замечание 1.1. В случае стационарного уравнения (т = 0) задача (1.1.6) принимает вид седловой системы (системы такого типа возникают и без использования смешанных постановок, например, в задаче Стокса для вязкой несжимаемой жидкости или в задачах фильтрации). Подробное изложение теоретических результатов для задач седлового типа содержится, например, в монографии Брецци-Фортена [46]. В данной работе будем предполагать, что коэффициенты достаточно гладкие и не будем останавливаться на условиях гладкости, которым должны удовлетворять коэффициенты уравнения для корректности постановки задачи.
Замечание 1.2. Необходимо отметить, что в отличие от классической постановки задачи, для смешанной постановки краевое условие Неймана для температуры (второе из условий (1.1.3)) становится главным, а краевое условие Дирихле (первое из условий (1.1.3)) - естественным.
Замечание 1.3. Мы сразу ограничились рассмотрением однородных краевых условий, так как наличие неоднородности в краевых условиях не приводит к существенному усложнению задачи. Учет неоднородных краевых условий Дирихле или Неймана сводится к добавлению дополнительного слагаемого в правой части (1.1.6). Также не представляет трудности рассмотрение третьего краевого условия (условие Робена, физически соответствующее конвективному теплообмену по закону Ньютона с окружающей средой на границе) вида аТ + bw n = д на части границы.
Замечание 1.4. Одним из основных преимуществ использования смешанной постановки в виде системы уравнений первого порядка (и использование смешанного МКЭ для аппроксимации по пространству) является то, что в таком виде непосредственно аппроксимируются законы сохранения. В данном случае - закон сохранения энергии (и закон Фурье, как определяющее отношение, связывающее температуру и тепловой поток). Построенные на основе такой постановки задачи численные методы будут полностью (глобально) консервативными, то есть в каждый момент времени будет выполнен сеточный аналог закона сохранения энергии. 1.2 Дискретизация по пространству
В данном пункте рассмотрим сразу трехмерный случай d = 3. Зададим в области Г2, вообще говоря, неравномерную параллелепипедальную сетку, состоящую из Nx Ny Nz ячеек z обозначим, соответственно, как hx = ХІ+\ — Xi,hJy = yj+\ — yj, и hkz = Zk+i — %k- Далее также будет часто использоваться максимальный шаг сетки h = maXij7k{h x,hJy,hk}. Введем естественную нумерацию ячеек сетки, а также их сторон, следующим образом. Ячейки сетки нумеруются в порядке увеличения индексов (x,y,z). Стороны ячеек, ортогональных направлениям ж, у и z - в порядке увеличения индексов х — у — z, y z xnz x y} соответственно.
Как уже упоминалось, для аппроксимации по пространству применяется смешанный метод конечных элементов, таким образом, неизвестные функции представляются с помощью МКЭ на каждом временном слое. Для аппроксимации температуры и других скалярных функций (плотности, теплопроводности и др.) будем использовать пространство кусочно-постоянных функций так что используемый метод конечных элементов относится к так называемым конформным методам конечных элементов [40]. Степени свободы для скалярных функций будем ассоциировать со значениями функций в центрах ячеек сетки, в то время как степени свободы для компонент векторных величин будем "привязывать" к серединам граней, ортогональных соответствующей оси координат. Таким образом, в терминах разностных схем мы используем аппроксимацию на разнесенных сетках. Замечание 1.5. На самом деле, в качестве степеней свободы можно использовать и интегральные средние, для скалярных величин - по ячейке, для компонент векторных величин - по соответствующим граням. Например, с точки зрения нахождения теплового потока в начальный момент времени как решения следующего уравнения в интегральной форме a(w,u) + 6(u,T0) = 0 VueV с заданной начальной температурой То при использовании элементов Ра-вьяра-Тома наименьшей степени требуются именно интегральные средние температуры по ячейке. Использование поузлового сброса значений в центры ячеек вносит дополнительную погрешность порядка /г2, что не влияет на порядок сходимости.
Полное изложение теоретических результатов, касающихся оценки погрешности интерполяции в нормах пространств Соболева для конечных элементов Равьяра-Тома можно найти, например, в [46], [57], [47]. Мы ограничимся здесь только повторением основного результата ([46], глава III):
Пусть Th - регулярное семейство триангуляции области Q и П - оператор глобальной интерполяции для элементов Равьяра-Тома наименьшей степени, Ну, V П (Ls(Q))d — V/j. Тогда существует константа с, не зависящая от h такая, что Кроме того, смешанный метод конечных элементов на параллелепи-педальных сетках обладает свойствами суперсходимости в сеточных нормах. Так, как следует из приведенных в [54], [55], [67], [63], [48] результатов, при наличии достаточной гладкости решения имеют место следующие оценки:
Схемы расщепления на основе скалярных схем расщепления для дивергенции теплового потока
В п. 2.1 представлены схемы расщепления с треугольной факторизацией оператора на верхнем слое по времени в уравнении для теплового потока. Данная идея рассматривалась в работах П.Н. Вабищевича [3], [4], [6] и в работах К.В. Воронина совместно с Ю.М. Лаевским [11], [14], [12]. В этом пункте рассматривается факторизация на основе попеременно-треугольного метода (пп. 2.1.1) и на основе метода симметричной верхней релаксации (SSOR) (пп. 2.1.2). Показана абсолютная устойчивость построенных схем для случая произвольной размерности.
В п. 2.2 рассматривается потоковая схема расщепления на основе алгоритма Удзавы [43] в двумерном случае. В п. 2.3 предложен новый подход к построению потоковых схем рас щеплення, основанный на использовании классических (скалярных) схем расщепления для сеточной дивергенции теплового потока. Предложенных подход позволяет построить большой класс схем, включащий в себя в качестве частных случаев схемы из п. 2.1-2.2. В подпункте 2.3.1 представлены потоковые схемы расщепления в двумерном случае, в подпункте 2.3.2 - в трехмерном случае. Общий результат, касающийся порядка аппроксимации потоковых схем расщепления, построенных с помощью предложенного подхода, сформулирован в пп. 2.3.3. В качестве порождающих схем для сеточной дивергенции используются схема переменных направлений, локально-одномерная схема, схема Дугласа-Ганна, а также схемы типа предиктор-корректор.
Наконец, пункт 2.4 главы 2 посвящен вопросу устойчивости потоковых схем расщепления на основе схем расщепления для сеточной дивергенции теплового потока. В пп. 2.4.1 сформулированы общие утверждения об устойчивости для температуры и устойчивости для теплового потока в полунорме, связанной с оператором дивергенции. Кроме того, в п. 2.4 для потоковых схем расщепления на основе схемы переменных направлений (пп. 2.4.2) и схемы Дугласа-Ганна (пп. 2.4.3) приведены полученные априорные оценки глобальной устойчивости по начальным данным и по правой части для температуры и теплового потока.
Основной идеей получения экономичных численных алгоритмов для решения рассматриваемой задачи является построение схем расщепления, так или иначе заменяющих обращение многомерного оператора G = A + 7jBM-1BT в уравнении на тепловой поток (1.4.8) на последовательное обращение некоторого числа одномерных. Среди способов построить схемы расщепления выделим следующие:
1. Приближенная факторизация оператора G = А + ВМ_1ВТ в явном виде. Например, при использовании попеременно-треугольной факторизации или факторизации типа SSOR оператор на верхнем временном слое заменяется на произведение блочно-треугольных операторов с трехдиагональными блоками на диагонали. Данный подход был рассмотрен П.Н. Вабищевичем и К.В. Ворониным совместно с Ю.М. Лаевским и описан в п. 2.1 главы 2.
2. Схемы на основе алгоритма Удзавы. Алгоритм Удзавы [53] можно понимать как обобщение метода переменных направлений на случай смешанной постановки. Данный подход был использован в работе Арбогаста и его соавторов [43] и описан в п. 2.2 главы 2.
3. Схемы расщепления, полученные с помощью предлагаемого в диссертации нового подхода, основаны на использовании классических (скалярных) схем расщепления для сеточной дивергенции теплового потока. Данный подход рассматривается в п. 2.3.
Основной идеей является приближенное представление оператора G на верхнем слое по времени в уравнении (2.0.1) в виде произведения блочно-треугольных множителей с блочной диагональю, которую можно обратить экономичным образом. При этом данный способ построения схем по времени может быть применен для произвольной размерности d = 2,3,..., поэтому мы не будем в данном пункте рассматривать отдельно
Априорные оценки и гладкость начальных данных
Заметим, что L = L(BTw) = BTJ-"(L)w. Более того, если J-(L\) и T{L i) определены корректно, то В2 J- {Li)J- {L2) = L\L i T. Таким образом, определив J {L) для всех операторов L в схеме для дивергенции, мы фиксируем, каким именно образом мы от схемы для дивергенции потока переходим к потоковой схеме. Примеры выполнения такого перехода для конкретных схем расщепления уже были представлены в пп. 2.3.1-2.3.2. Например, легко проверить, что действие J- на основные операторы (напомним, в пункте 2.3 были введены операторы Ах = В ВХ и Ах =
Замечание 2.5. (о "подходящих" операторах L) Так как на части границы ставится условие Дирихле, то либо оператор Лж, либо Ау является невырожденным. Поэтому, если, например, Ах - невырожденный, то можно записать произвольный оператор L в виде L = АХАХ L , формально выполнить описанную выше процедуру и определить
Отсюда следует вывод, что формально оператор L может быть произвольным и что, используя разные представления одного и того же оператора, можно по-разному определять J (L): то есть J- не является однозначным отображением. На самом деле, если в оператор L\ (или L2,iv3,iv4) входят обратные операторы Л"1 или Л"1, то обращение L\ (или L2,iv3,iv4) уже не может быть выполнено экономичным образом, так как обратные операторы требуют обращения заполненных матриц (вида ВХА ВХ} где Вх и Ах имеют вид 1.3.4 или 1.3.5 и 1.3.6 или 1.3.7). Напомним, локально одномерная схема для сеточной дивергенции теплового потока может быть записана в следующем виде (2.3.4): которая имеет в перестановочном случае (ЛжЛу = ЛуЛж) второй порядок точности, в неперестановочном - только первый. Соответствующая схема расщепления для теплового потока w тогда выглядит (с учетом упрощенных обозначений (1.4.11)) следующим образом: так что равенство невозможно, и в схеме для потока соответствующие операторы неперестановочны. Из этого следует, что порядок аппроксимации у потоковой схемы расщепления на основе локально одномерной схемы для дивергенции будет таким же, как и у локально одномерной схемы для дивергенции в неперестановочном случае (то есть - первым). Этот факт можно обобщить до следующего утверждения: Порядок аппроксимации потоковой схемы расщепления не ниже порядка аппроксимации порождающей схемы расщепления для дивергенции теплового потока в неперестановочном случае.
В данном пункте будут рассмотрены теоретические результаты по устойчивости для некоторых потоковых схем расщепления, построенных с помощью разработанного в пункте 2.3 подхода для двумерного и тремерно-го случаев. Общие результаты по устойчивости потоковых схем расщепления приведены в п. 2.4.1. А именно, получены оценки устойчивости по начальным данным и по правой части для температуры, а также устойчивость в полунорме для теплового потока. В пи. 2.4.2 будут приведены априорные оценки устойчивости для потоковой схемы расщепления из п.2.3.1 на основе схемы переменных направлений для сеточной дивергенции теплового потока (двумерный случай), а в подпункте 2.4.3 - для потоковой схемы расщепления из п.2.3.2 на основе схемы Дугласа-Гана (трехмерный случай). Всюду в данном пункте используются упрощенная форма записи (1.4.11). 2.4.1 Устойчивость в полунорме для потоковых схем расщепления. Общие результаты
Напомним, в (1.4.13) в пункте 1.3 главы 1 были введены конечномерные гильбертовы пространства Н и Щ (сеточные аналоги пространств V и Q7 см. (1.1.5)) для теплового потока и температуры, соответственно. Основная идея получения априорных оценок состоит в использовании возможности перехода от схемы в целых шагах для дивергенции теплового потока к схеме в целых шагах для теплового потока. Пространство Н представим в виде прямой суммы сеточно-соленоидальных и сеточно-потенциальных векторов: эвклидова норма в H. В соответствии с подходом из пункта 2.3, предположим, что потоковая схема расщепления была получена из некоторой скалярной схемы расщепления для сеточной дивергенции теплового потока = BTw. Предположим также (для каждой конкретной порождающей схемы расщепления для дивергенции данное условие можно проверить), что схема для устойчива в сеточном пространстве с нормой II и, более того, имеет место следующее неравенство:
Например, в перестановочном случае обычно имеет место устойчивость схемы расщепления для в сеточном пространстве, эквивалентном пространству L И И = 2 В соответствии с (1.4.7) и, учитывая упрощенные обозначения (1.4.11), при / = 0 имеем:
Тогда в уравнении для теплового потока (например, для схемы Кран ка-Николсон см. (1.4.6)) появится правая часть Fn, где Fn = Brn Є /mB. ЭТО означает, что если от уравнения в целых шагах для можно перейти к уравнению в целых шагах для соответствующей потоковой схемы расщепления, то в уравнении для сеточной дивергенции правая часть будет иметь вид: sn = Arn . Таким образом, вопрос об устойчивости по правой части для температуры сводится к известным результатам для схемы расщепления для сеточной дивергенции теплового потока.
Замечание 2.6. В данном пункте были приведены только соображения, приводящие к устойчивости теплового потока в полунорме, связанной с оператором дивергенции В . Для потоковых схем расщепления в двумерном и трехмерном случаях, основанных на схемах переменных направлений и схеме Дугласа-Ганна, априорные оценки глобальной устойчивости по начальным данным для теплового потока будут приведены в двух следующих подпунктах.
Моделирование термохронологии косой позднепалеозойской коллизии в зоне Таймыр Северноземельской складчатости
В этом пункте рассматриваются термохронологичекие реконструкции на этапах магматической эволюции в неопротерозойской активной континентальной окраине Сибирского кратона. Основной целью было численное моделирование термальных режимов в рассматриваемой области и проверка ряда геологических гипотез [78]. Пп. 4.1.1 и 4.1.2 содержат описание физической и математической моделей соответственно. Задача рассматривалась в двумерной постановке, для численного моделирования применялась двумерная потоковая схема расщепления на основе схемы переменных направлений для сеточной дивергенции теплового потока (см. п. 2.3). В пп. 4.1.3 приведены результаты численных экспериментов и кратко сформулированы основные выводы, полученные в ходе решения данной прикладной задачи.
На Рис. 4.1 изображена рассматриваемая физическая область (вид сверху), отрезок АА обозначает линию пересечения с плоскостью 2D-MO-делирования. Ориентировка модельной области соответствует разрезу западно-восточного простирания через область Татарско-Ишимбинской шовной зоны. После существенного упрощения реальных геологических данных и перехода к двумерной постановке, модельная область выглядела так, как показано на Рис. 4.2.
Модельная область (схематично) с указанием основных пород. 121 км - в ширину. На нижней границе области был задан однородный тепловой поток wn = — 0.05кг , направленный внутрь области, на верхней границе - заданная температура 300К (температура на поверхности Земли). Так как боковые границы области появились в результате искусственного ограничения рассматриваемой области, для простоты на этих границах было поставлено условие равенства нулю теплового потока через границу. При этом слева и справа в область модели были добавлены "буферные зоны" с целью уменьшить влияние краевых условий на распределение температуры в центральной части модели.
Вмещающие породы образовывали периодическую структуру и состояли из слоев со слюдисто-карбонатным сланцем, кварц-биотитовым сланцем, кварцитом и кальцитовым мрамором (см. подписи на Рис. 4.2). В начальный момент времени в области также находились амфиболитовые тела (обозначены зеленым цветом). Кроме того, в соответствии с описанным ниже временным сценарием, в области модели в определенные моменты времени "мгновенно" появлялись тела карбонатит-гранитной ассоциации - щелочное габбро и карбонатиты пенченгинского комплекса (вытянутые тонкие тела оранжевого и бледно-розового цветов), а также лейкограниты (красный) и граниты (ярко розовый) Татарского массива. В таблице 4.1 приведены физические характеристики рассматриваемых типов пород.
Далее приведено краткое описание основных этапов моделирования. Рассматривался промежуток времени длиной 125 млн. лет, начальный момент времени - 725 млн лет. назад. В определенные моменты времени, согласно геологическим данным, в рассматриваемую область внедрялись тела карбонатит-гранитной ассоциации. При этом с точки зрения геоло 122 гических времен это внедрения происходило очень быстро, практически "мгновенно". В связи с этим процесс моделирования был разделен на несколько этапов, связанных с внедрением и последующим остыванием различных тел. В начальный момент времени (725 млн. лет назад) произошло внедрение карбонатитов и габбро, далее через 79 млн. лет (646 млн. лет назад) в области модели появились лейкограниты, и, наконец, еще через 17 млн. лет произошло внедрение гранитов (629 млн. лет). Очевидно, основные изменения температуры в области были связаны с наиболее объемными интрузивными массивами гранитов и лейкограни-тов.
Одной из особенностей расссматриваемой модели являлся учет флю-идонасыщенности на этапах внедрения гранитов и лейкогранитов. Влияние этого фактора на рассматриваемую тепловую задачу сводилось к тому, что в некоторой зоне вокруг гранитов и лейкогранитов был повышенный тепловой поток (за счет введения тензора эффективной теплопроводности) вокруг внедренных массивов. Как показывают результаты, приведенные в пи. 4.1.3, наличие флюидонасыщенных пород оказывает заметное влияние на распределение температуры в центральной части модели.