Параллельные технологии решения краевых задач

Василевский, Юрий Викторович

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

Василевский, Юрий Викторович АВТОР

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

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

2005 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

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

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

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

Василевский Юрий Викторович

ПАРАЛЛЕЛЬНЫЕ ТЕХНОЛОГИИ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ

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

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

Москва - 2006

Работа выполнена в Институте вычислительной математики РАН

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

профессор Агошков В.И.

доктор физико-математических наук, профессор Лапин A.B.

доктор физико-математических наук, профессор Чижонков Е.В.

Ведущая организация: Институт математического моделирования

РАН

Защита состоится 27 апреля 2006 г. в 11 часов на заседании диссертационного совета Д 002.45.01 в Институте вычислительной математики РАН по адресу: 119991, г. Москва, ГСП-1, ул. Губкина, 8.

С диссертацией можно ознакомиться в библиотеке Института вычислительной математики РАН.

Автореферат разослан 24 марта 2006 г.

Ученый секретарь диссертационного совета доктор физико-математических наук

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

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

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

Цель исследования заключается в конструировании и исследовании новых параллельных численных методов адаптивного решения краевых задач с использованием разных типов сеток (анизотропных неструктурированных, нестыкующихся, иерархических) и применении этих методов к решению актуальных прикладных задач.

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

Научная новизна и теоретическая значимость. В рамках теории конечно-элементных

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

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

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

Апробация работы. Основные результаты диссертации докладывались автором: на всероссийской школе-конференции (Казань, 1999), всероссийских конференциях по построению расчетных сеток (Москва, 2002,2004), международных конференциях по методам декомпозиции области DDM (Чиба, 1999), параллельной вычислительной гидродинамике ParCFD (Москва, 2003), геофизическим наукам SIAM (Остин, 2003, Авиньон, 2005), Европейских конференциях по вычислительной математике ENUMATH (Париж, 1995, Юваскюла, 1999, Иския, 2001), Европейской конференции по вычислительным методам ECCOMAS (Париж, 1996), Франко-Русско-Итало-Узбекских симпозиумах по численному анализу и приложениям (Ташкент, 1995, Марсель, 1997), международной конференции по анализу, вычислениям и применениям дифференциальных и интегральных уравнений (Штуттгарт, 1996), международной конференции GAMM по параллельным многосеточ-'ным методам (С.-Вольфганг, 1996), международной конференции "Domain decomposition and multifield theory" (Обервольфах, 1998), международных симпозиумах "Finite element workshop" (Хьюстон, 1999, Колледж Стэйшн, 2002) международной конференции "50 лет сопряженным градиентам" (Юваскюла, 2002), на научно-исследовательских семинарах Института вычислительной математики РАН, Вычислительного центра РАН, Института математического моделирования РАН, Университетов Париж 6, Париж 13, Лион 1, Ренн, Гейдельберг, Мюнхен, Аугсбург, Юваскюла, Наймеген, Остин, Хьюстон, Национальной Лаборатории в Лос Аламосе, а также на семинарах INRIA, Institut Francais du Petrol, ExxonMobil Upstream Research C.

Публикации. По теме диссертации опубликовано 42 работы, из них 21 в рецензируемых журналах, 18 в материалах конференций, 3 в научных изданиях.

Личный вклад автора. Вклад автора в совместные работы заключался: в формировании постановки проблемы [2,3,4,5,6,7,8,9,10,11,12,16,17,18,19,20,21|, предложении идеи решения [2,3,4,5,6,7,8,9,10,11,12,13,14,16,17,18,19,20,21],теоретическом обосновании [3,4,6,10, 11,16,17,20], совместном теоретическом обосновании [2,5,7,8,9,12,13,14,15,19,21], технической реализации [6,13,14], совместной технической реализации [3,5,7,9,10,11,12,15,16,17,18, 19,20,21], постановке численных экспериментов [3,5,6,7,9,10,11,12,13,14,15,16,17,18,19,20,21]

Структура диссертации. Диссертация состоит из введения, четырех глав и заключения. Текст работы изложен на 254 страницах, содержит библиографию из 275 наименований, 33 рисунка и 57 таблиц.

Содержание диссертации

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

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

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

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

Далее рассматри ваются новые комбинации известных методов, оказывающиеся весьма эффективными для неструктурированных систем специального вида. Быстрое приближенное решение систем с матрицами, порождаемыми аппроксимациями на неструктурированных призматических сетках, рассмотрено в разделе 1.2. Такие матрицы обладают тензорным рангом 2, т.е. представимы в виде >11 ® Мчз + М\ ® Лгз- Здесь А?з, Л/23 М1) — матрицы жесткости и масс (диагонализированная) на двумерном (со-

отв. ,одномерном) сечении призматической сетки. Метод представляет собой обобщение быстрого прямого метода (Ю.А.Кузнецов) на случай неструктурированных матриц Л23 с двумерным графом и использует явный вид матриц А\, Лгз, Мц- Смысл обобщения — замена прямого решения двумерных систем с матрицами Л23 + ХМ23 с произвольным параметром Л на приближенное итерационное решение с переобуславливателем, базирующемся на факторизации одной из трех матриц Л23+ЮМ23, Л23 + 150Мгз, А23 + 15ООМ23. При этом разрешающий оператор будет нелинейным, и для решения системы применим быстро сходящийся итерационный метод FGMR.ES с построевшым нелинейным переобуславливателем, арифметическая цена применения которого на практике пропорциональна порядку матрицы с полилогарифмическим множителем.

В разделе 1.3 излагается новая технология переупорядочивания неизвестных, позволяющая для диффузионных задач с анизотропными коэффициентами успешно использовать блочные метод Гаусса-Зейделя и легко параллелизуемый метод Якоби. Даже для аппроксимаций на неструктурированных сетках удается переупорядочить матрицу и разбить ее на блоки таким образом, что эти методы становятся весьма эффективными. Для

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

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

Рассмотрим конформную триангуляцию Пл области Г2, которая разбита на т перерывающихся подобластей регулярной формы с минимальным перекрытием в один элемент. Пусть матрица А получена стандартной конечно-элементной Р\ -аппроксимацией диффузионного оператора с коэффициентом р(х) с условиями Дирихле. Разбиения сетки и области порождают блочное представление матрицы А, для которого можно определить блочно-диагональный переобусдавливатель В\:

Здесь блоки Вц = В? > 0 размерности П; — переобуславливатели для блоков Ац, г = 1,..., т, с константами эквивалентности 0 < ац < а?. Матрица — легко парал-лелизуемый аддитивный переобусдавливатель Шварца с минимальными перекрытиями между подобластями. Его эффективность зависит от числа подобластей т, ширины налегания 5 и скачка диффузионного коэффициента р. Для того, чтобы уменьшить негативное влияние малого налегания и исключить зависимость от числа подобластей и скачков р(х), применим поправку на некотором грубом сеточном уровне. Результирующий гибридный декомпозиционный переобусдавливатель Вы определяется выражением

где матрица Вц задается следующим образом. Пусть число ненулевых строк в матрице А = [-<4<1, ■ ■ •, А,«+1) •■• 1 Д,г»] (без диагонального блока ) равно гч. Определим п =

Щ + т и предположим, что строки матрицы А упорядочены таким образом, что в

каждой матрице Л; ненулевые строки идут первыми. Тогда локальная Тц и глобальная

Ви = (/ - В2А)В1(1 - АВ2) + В2,

Т матрицы агрегирования задаются как

Ти= f J ° I , Г= | ••. | , Тц € Rn'x(ii|+1)

где = (1,..., I)7" 6 R"1-"', — единичная матрица.

Определим на "грубом" подпространстве агрегированную матрицу жесткости А — ТтАТ 6R"™ и переобуславливатсль В для нее. Матрица Bi = ТВТ1" — переобуслав-ливатель в пространстве агрегированных векторов К = R": v = Tv, v £ R"}.

Теорема 1.3. Пусть f2 разбита на m налегающих подобластей регулярной формы Çli с шириной минимального налегания S и диаметром Н, триангуляция П квазиравномерна в области налегания, расширенной на один сеточный слой, коэффициент р(х) гетерогенен в области налегания подобластей и гладок во внутренних (неполегающих) частях подобластей: рц, < р(х) < p2,i, х € и пусть В и Вц, г = 1,... , m, — переобуславливатели для А и Ац, соответственно, с константами эквивалентности О < c*i < ос2. Тогда для числа обусловленности матрицы В/,А верна оценка:

Concï(BhA) < 1— (l + 2С тах (• - ai V VPi.i/ <V

Здесь и далее с, С (с индексами и без индексов) будут обозначать некоторые положительные константы, не зависящие от параметров задачи, а выражение А ~ В будет обозначать спектральную эквивалентность матриц Л и В с константами, не зависящими от размеров матриц, или приближенное равенство или пропорциональность величин А и В.

На практике оказалось, что число обусловленности не зависит от гетерогенности коэффициента не только в зоне налегания, но и внутри подобластей. Отметим, что можно использовать более экономичный вариант переобуславливателя В = В\ 4- Вг(/ — АВ1), который при выборе В = А'1 и соответствующей коррекции начального приближения, щ := Ио + Вг(/ — Лио), оказывается эквивалентным В*,.

Тип переобуславливателя В зависит от конкретного приложения. В двумерном случае достаточно использовать прямую факторизацию матрицы Л, поскольку ее размерность не велика, а методы факторизации являются весьма эффективными для матриц с двумерными графами, как показано в разделе 1.1. В трехмерном случае можно использовать методы неполной факторизации. При параллельных расчетах для переобуслав-ливапия А удобно применить несколько итераций ВБОК, где блоки получены жадным алгоритмом расцвечивания и переупорядочивания строк и столбцов. Последний подход

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

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

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

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

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

Определение 2.1. Пусть задана и 6 С°(П) и оператор С°(П) —» Pi(Clh) — проек-

тор па пространство непрерывных кусочно-линейных функций .Pi(Па) • Сетка fiJ^/Vr, и), состоящая из не более чем Nt элементов, называется оптимальной, если она — решение оптимизационной задачи

и) = arg ш llu-'PsJ^jU^n).

Здесь J\f(flh) обозначает число симплексов в Пл.

В качестве оператора Vnh могут выступать, например, оператор Pi -интерполяции Хп„ или конечно-элементная аппроксимация краевой задачи.

Большинство теоретических результатов, представленных ниже, основано на предположении, что функция и 6 C2(fi), a Vnh — ~ оператор интерполяции. Однако постоянные, входящие в наши оценки ошибок, не зависят от действительного значения С2-нормы и. Это позволяет применять оценки для функций, близких к функциям с характерными для краевых задач особенностями. Если Vah ф Та^. то теоретическое обоснование базируется на связи между ||u — Pn^Hi»^) и ||и — Хг^иЦ^п).

При определенных предположениях на оператор VnЛ можно доказать существование оптимальных сеток в двумерном случае.

Теорема 2.3. Пусть и € С°(П) и ||ы — 7-Vikullioo(n) — непрерывный функционал от координат сеточных узлов, т.е.

III« - VnM\Loc(ii) - II" - ^«11^(0)1 < C(u)£,

где HJ — триангуляция, полученная любым £ -возмущением узлов некоторой конформной триан.гуляции П),. Кроме этого, пусть проектор Vak удовлетворяет неравенству

II« - ■Pnjullwri) 2 II« - ^n'«IU„(n)

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

Если существование оптимальной сетки доказано только в двумерном случае, то ее свойства проанализированы как в двумерном (<i = 2), так и в трехмерном (d = 3) случаях.

Теорема 2.4- Пусть гессиан Н функции и € С2(П) невырожден в П, и пусть для любого симплекса А 6 П°,'* выполняется следующая оценка:

И-Яр, ~ (ЯдЫ^д) < 9д1МЯд)|, О <<?д < q < 1, p,s =

где Яд = H(arg max|(2etH(z)|), Ах(Яд) — ближайшее к 0 собственное значение Яд. Тогда

Здесь |П||я| обозначает объем области П в метрике |Я|, описанной ниже.

Верхняя граница для ||и — Гп^мЦ^р) аналогична нижней и доказана в Теореме 2.8 через оценку ошибки на квази-оптимальных сетках. Определение оптимальных сеток допускает обобщение на случай Lv-нормы:

Определение 2.6. Пусть заданы и € С°(П) и р

€]0,+оо]. Сетка U^iNr,

и), состоящая из не более, чем Nt элементов, называется оптимальной по отношению к Lp-чорме, если она — решение оптимизационной задачи

ПT(NT, и) = arg min ||u - *>Пк«||Мп,.

Теорема 2.7. Пусть выполнены условия Теоремы 2.4- Для случая незнакоопределенно-го гессиана Н предположим дополнительно существование положительных констант со,с\, таких, что для любого симплекса ■ Д С П^1 верно

С1(ЯдС,С) УС 6 НЛ Тогда для оптимальных по отношению к Ьр-норме сеток справедлива оценка

С(Я, со)^^ < II" - ,Я| := (сге(|Я|)-1/(2р+<')|Я|.

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

ности в заданной непрерывной метрике С. Пусть |е|с — объем симплекса е в метрике £7, а \8де\а — суммарная длина его ребер в метрике С. Тогда 0 < < 1 задается

через качество элементов <3с(е):

^(ад-шцдоМ. <?.(*) =

где

и |Г2д|<з = 1е|с — объем расчетной области в метрике <7. Таким образом, если качество

сетки достигает 1, то она состоит из равносторонних в метрике С симплексов диаметра Л*. Сетки с качеством <?с(Пл) ~ 1 называются (З-квази-равномерными.

Пусть Я — невырожденный гессиан фунцкции и б С2(П) и Я = — его

спектральное разложение с ортонормальной матрицей IV и диагональной матрицей Л = Лар{Л;}, г = 1,...,с2, (|Л1| < ••■ < Определим метрический тензор (в дальней-

шем просто метрика) |Я| = "[^[^У/. Основным конструктивным результатом является утверждение о том, что |Н| -квази-равномерная сетка является квази-оптимальной, т.е. аппроксимирует оптимальную сетку.

Теорема 2.8. Пусть гессиан Я функции и € С2(П) невырожден е П, и пусть — \Н\-квази-равномерная сетка с /Уг элементами, такая, что £?|#|(^л) > (¿о- Далее, пусть для симплекса Д 6 Пл, где достигается ||и — 2ЬЛи||£„,(а), выполняется оценка:

||яра - (яд)р„||^(д) < <7д|м#л)1. о < Ял < я < 1, р,5 =

где Яд = Я(aгg тад \<1е1Н{х)\), А^Яд) — ближайшее к 0 собственное значение Яд. Тогда П/, является квази-оптимальной сеткой, т.е.

11« - 2п„и|к=о(п> £ СЮо,?)!^ -Х^р.иЦ^да,

причем ||и-Хп„Чи.<п) <С(С?о,<?) (^г) ' ■

Для проекторов ~РПн, удовлетворяющих ||и - Р^иЦ^^) < С|[и - ХиЦ^п), оценка ошибки преобразуется к виду

И»-■ад< сс(д0,9)

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

Как правило, гессиан Н(х) является неизвестной тензорной функцией. На практике используется его приближение Я|> = {Я£,}р<=1) Яр, € РДП/,), восстановленное в узлах сетки из сеточной функции ик = Тпки € Р^Пл) следующим образом. Для внутреннего узла а;

У Я* (а>к Ах = - I ^ ^ 6 Р1(<т0. ^ = 0 "а 3<т;,

СГ»

где суперэлемент cr¡ — объединение симплексов с общим узлом щ. В граничных узлах Hp,(a¡) — взвешенные экстраполяции значений с близлежащих внутренних узлов.

Теорема 2.11. Пусть и £ CJ(f¡), uh = Vi\hu, Н — гессиан и, и Hh — восстановленный дискретный гессиан. Более того, пусть для любого суперэлемента а € Пл, ассоциированного с каким-то узлом а, выполняются следующие оценки:

\\НР, - На,р,< 6, |tf¿» - H„,paI < e,

где Ha = H(x„) и х„ = arg max |det(ff(x))|). Тогда для е и 5 достаточно малых. по отношению к минимальному собственному значению \На\, \Нк\-квази-равномерная сетка П/, ä Qo) является также и \Н\-квази-равномерной: >

С Qo где постоянная С не зависит от Nj- и и.

Теорема 2.11 утверждает, что при некоторых предположениях достаточным условием квази-оптимальности является |./У'1|-квази-равномерность. Первое предположение означает малые вариации гессиана на любом суперэлементе сг, что достигается на сетках с большим JVr, а второе предположение — требование аппроксимации гессиана в узлах, что подразумевает малую градиентную ошибку для ин. Малая градиентная ошибка не характерна для решений с особенностями. Для восстановления дискретного гессиана в случае негладких функций мы предлагаем другое определение его компонент удовлетворяющих второму предположению в более слабой норме.

Определение 2.12. Триангуляция удовлетворяет условию А, если для любого внутреннего суперэлемента существует аффинное отображение .F; = <S¿ ° такое, что — регулярный суперэлемент диаметра 1, для которого радиус наибольшей вписанной сферы 0(1). Здесь St и обозначают матрицы масштабирования и вращения, соответственно.

Фактически, условие А — требование локального подобия форм соседних симплексов. Пусть для узла at сетки, удовлетворяющей условию A, обозначает наибольшую сферу с центром J~i(ai), вписанную в В силу регулярности радиус R¡ сферы

0(1). Вводя полярные координаты с центром в Pifa), мы можем определить гладкую функцию V; = 1 —т2/Щ на B¡. Линейная оболочка функций v = aj7^1^), а 6 R1, задает пространство пробных функций V¡. Отметим, что v € Vi означает v 6 C2(B¿), v = 0 на dB¡, где B¡ = |В<| > |ст;|. Для повернутого симплекса Дтг = TIA

можно определить его размеры h-^ f¡ = max |(i)t — (2/)¡t|, k = 1,..., d.

Во внутреннем узле сетки a¡ компоненты удовлетворяют равенству

/ = / - J иНкпЖ VlJ"€

B¡ В, OBi

В граничных узлах H^Jfli) ~~ взвешенные экстраполяции значений с близлежащих внутренних узлов.

Теорема S.J5. Пусть заданы функция и € П IV^fi), р > d, и внутренний

суперэлемент , удовлетворяющий условию А . Предположим, что на симплексах Л; С ст, градиентная ошибка для функции наилучшей линейной аппроксимации Як функции u-ii — изотропно распределена:

majc/i£||0"(u*. -Ш)\\ьР(дя) < Св min Л|||0°(иге - E£)||t

|о| = 1 М—1

причем имеет место сходимость «к к :

h°n\\d"{un - й£)||Мдя) < Cch^a0, |а| — 1, 0 > О,

где 0 — параметр дополнительной гладкости и. Кроме того, пусть на <Tj дифференциальный гессиан Н отклоняется незначительно от своего среднего Н: ||НГа—¿¿»lUito) £ <5, г, s = 1,..., d. Тогда дискретный гессиан НТ,, восстановленный из кусочно-линейного интерполянта сходится к дифференциальному гессиану:

II Hr. - нг\ ||WBi) < S + CbCcW^

причем оценку можно упростить: ||НТ, — ~ ^ + CßCc

Отметим, что последняя оценка означает локальную сходимость дискретного гессиана в слабой норме, т.к. для любой триангуляции, адаптированной к функции с невырожденным гессианом, lim max min hv k = 0- Важность теоремы в том, что это -v(su)-. те k=l,...,d

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

Пусть G — непрерывная кусочно-линейная метрика, заданная через значения своих компонент в узлах сетки. Построение G-квази-равномерной триангуляции, состоящей из приблизительно Nt симплексов, — ключевой элемент технологии адаптивного построения сеток на основе восполнения дискретного гессиана сеточного решения. Симплексы сетки должны как можно меньше отличаться от равностороннего (в метрике G ) симплекса с длиной ребра h*, чей объем равен \CI\o/Nt- Алгоритм состоит в генерации

г I

последовательности сеток flj"™,..., П^",..., П\, такой, что

Qc.Nri^l) < Qo.nAnütö < • < Qg.nA^) <■■■< Qg.nASÜ).

Каждая сетка — член последовательности — является локальной модификацией предыдущей сетки, не уменьшающей ее качества по отношению к h* в G и представляющей

собой следующие действия. Сначала выбирается симплекс с наихудшим качеством. Затем формируется множество симплексов текущей сетки, имеющих ненулевое пересечение с выбранным симплексом. К созданному множеству пробуется ряд топологических операций, которые сохраняют его границу и могут повысить локальное качество сетки. Та виртуальная операция, которая повышает локальное качество, выполняется реально и приводит к очередной сетке ПД"'" . Если ни одна из операций не повышает качества сетки, то рассматриваемый симплекс на некоторое время удаляется (фиктивно) из сетки, и операции модификации пробуются к множеству, ассоциированному со следующим симплексом наихудшего качества. Подобная гибкость обеспечивает устойчивость алгоритма к появлению неразрешимых ситуаций. Важным технологическим преимуществом этого алгоритма является то, что он может быть реализован как черный ящик.

Алгоритм генерации квази-оптимальных сеток

Шаг инициализации. Построить начальную сетку О.^. Выбрать качество конечной сетки <5о, <2о < 1, и нужное число УУг сеточных элементов.

Итерационный Шаг.

1. Вычислить сеточное решение ил = "Рпк«•

2. Восстановить дискретный гессиан Нк сеточной функции и*1. Если > <Эо, остановиться.

3. Построить следующую сетку ¿1/, такую, что

4. Положить П/1 := П/, и перейти к 1.

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

Вопрос анализа сходимости адаптивного алгоритма для решения краевых задач с особенностями остается открытым. Для случая минимизации ошибки кусочно-линейной интерполяции гладкой функции рассмотрен механизм сходимости адаптивного алгоритма. Теорема 2.18 утверждает, что при определенных условиях в ходе адаптивных итераций градиентная ошибка уменьшается и восстановленный гессиан сходится поточечно к дифференциальному, обеспечивая рост качества сетки Я\ц\,нт{&и) и, следовательно, уменьшение интерполяционной ошибки.

Одно из основных преимуществ представленного алгоритма — простота управления адаптацией за счет модификации метрики — рассмотрено в разделе 2.3. Исследуются модификации \Н\ метрики |Я|, основанные на модификации собственных значений в ее спектральном разложении: |Я| = WT\K\W, \Н\ = WT\A\W.

Теорема 2.19 Пусть и € С2(П), гессианы Н и Н — невырождены в П, Г2ь — |Я'Ч -квази-равномерная сетка, Q|/7h|(S^) > Qo, а П/, — \Hh\-квази-равномерная сетка, <3|яч(Г2л) > Qo- Кроме того, пусть для любого симплекса Д € ÍÍ/, и симплекса Л* 6 П/t, где достигается ||u — Xñ^ilU^tn) > выполнены следующие оценки:

||ЯР, - (ЯЛ)р,|и„(д) < ддА^ЯдШ, ||Ярз - (Яд)р,||^(Л) < íaA^I^DI,

для 0 < <?л < <7 < 1 и р, 8 = Если Я£. — незнакоопределенный тензор, потре-

буем дополнительно, чтобы Сопс1(Лд-Лд1) q¿. < q. Тогда

< a(H,H)\\u-InKu]\L„{nh а(Н,Н) = C(Q0lq) mMc^lAA"1!) (pjj) '

где р(А) обозначает спектральный радиус матрицы А.

В качестве примеров метода управления рассматривается построение изотропной метрики, для которой A i = argmax{|Aj|,..., |A¿|} и Сопс1(|Лд.Лд!|) = j Ах |, а также ослабление (усиление) влияния особенности за счет весовой функции ш(х) > 0, при котором А< = Ajcj(i), i = 1,... ,d, и Со1к1(|Лд-Лд!|) = 1. В последнем случае верно Следствие 2.20. Пусть w < 1. Тогда верна следующая оценка:

II"-^"llí.S.OT ^ c(Qo,9)||"-Tf¡,,u||l,„,(n), ll"IUs,(n) :=sup|w(x)u(i)|.

ien

Таким образом, локальное подавление гессиана и результирующая релаксация сгущения адаптивной сетки приводят к оптимальным оценкам в весовой норме || •

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

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

собой следующие действия. Сначала выбирается симплекс с наихудшим качеством. Затем формируется множество симплексов текущей сетки, имеющих ненулевое пересечение с выбранным симплексом. К созданному множеству пробуется ряд топологических операций, которые сохраняют его границу и могут повысить локальное качество сетки. Та виртуальная операция, которая повышает локальное качество, выполняется реально и приводит к очередной сетке П,"'" . Если ни одна из операций не повышает качества сетки, то рассматриваемый симплекс на некоторое время удаляется (фиктивно) из сетки, и операции модификации пробуются к множеству, ассоциированному со следующим симплексом наихудшего качества. Подобная гибкость обеспечивает устойчивость алгоритма к появлению неразрешимых ситуаций. Важным технологическим преимуществом этого алгоритма является то, что он может быть реализован как черный ящик.

Алгоритм генерации квази-оптимальных сеток

Шаг инициализации. Построить начальную сетку Пд. Выбрать качество конечной сетки Цо, <Зо < 1, и нужное число N1 сеточных элементов.

Итерационный Шаг.

1. Вычислить сеточное решение ил = Тпьи.

2. Восстановить дискретный гессиан Нк сеточной функции и1*. Если <3|ял||^1.(Г2л) > <2о> остановиться.

3. Построить следующую сетку П/, такую, что

Я\н*\,ыт№ь) > За-

4. Положить П^ := Г2Л и перейти к 1.

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

Одно из основных преимуществ представленного алгоритма — простота управления адаптацией за счет модификации метрики — рассмотрено в разделе 2.3. Исследуются модификации |Я| метрики |Я|, основанные на модификации собственных значений в ее спектральном разложении: |Я| = WT\K\W, |Я| = WT\K\W.

Теорема 2.10 Пусть и 6 С2(П), гессианы Н и Н — невырождены в П, Clf, — \Hh\-квази-равномерная сетка, > Qo> а ilh — \Ны\-квази-равномерная сет-

ка, (5|я1>|(Г2л) > Q0- Кроме того, пусть для любого симплекса Д € Г2/, и симплекса А* ё ii/n где достигается ||и — Xft,,uIU«>(i4, выполнены следующие оценки:

||я„, - (Яд^Ц^^) < -гдЛ.аЯд!)!, ||ЯР> - (Яд)р,|и„(д) < ддА^ЯдШ,

для 0<дд<д<1 и р, s = 1,... ,d. Если Я£. — незнакоопределенный тензор, потребуем дополнительно, чтобы СопсЦЛд.Лд!) ?д- < q. Тогда

< a(H,H)\\u-lnhu\\L„in), a{H,H) = C{Q0,q) ткр(|ЛА-»|) (щ^) ' ■

где p{A) обозначает спектщльный jmduyc матрицы А.

В качестве примеров метода управления рассматривается построение изотропной метрики, для которой А; = argmax{|Ai|,..., |А^|} и Соп<1(|Лд-Лд1|) = |Arf/Ai|, а также ослабление (усиление) влияния особенности за счет весовой функции ш(х) > 0, при котором Ai = i = 1,... и Со1гс1(|Лд-Лд1|) = 1. В последнем случае верно

Следствие 2.20. Пусть ш < 1. Тогда верна следующая оценка:

IIй-:z'i!fcullx.Si.(n) ^ С(<2о,7)||и-ХП).и|и.»(п), IMUütn) :=sup|w(i)u(a;)|.

хеГ2

Таким образом, локальное подавление гессиана и результирующая релаксация сгущения адаптивной сетки приводят к оптимальным оценкам в весовой норме || • Hi^n).

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

Квадратичная экстраполяция Iрг функции ц>к с гессианом Я4*2 использует многоточеч-

ную формулу Тейлора = ХХЯ^-сц), (?-а()) где аь а2, а3 —вершины

, а — кусочно^линейная функция, такая что = <5^.

Теорема 2.22. Пусть ^ — квази-равпомерная триангуляция с тагом к, Я* и ЯЛ — дифференциальный и дискретный гессианы функций <р и цз^ соответственно, и

\\Щ. - < <5. И^(¥> - <

Тогда квадратичное восполнение функции <рь,, удовлетворяет оценке

где константа С не зависит от 5, е, к и у.

В частности, если <р принадлежит С3(ст'), то е ~ Л1, б ~ Л и оценка аппроксимации границы становится — ^гН^^гч) — что на порядок выше, чем у ¡рн-— ^аII¿^(г,) — С^2- В силу этого, влияние неточности задания границы на аппрок-симационную ошибку существенно уменьшается.

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

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

Использование неконформных аппроксимаций дает дополнительные возможности для построения параллельных адаптивных технологий решения краевых задач. В третьей главе, состоящей из пяти разделов, рассматриваются блочные параллельные методы решения систем линейных уравнений, возникающих в некоторых прикладных неконформ-иых аппроксимациях, и основанные на этих методах параллельные адаптивные технологии. Разбиение области на подобласти и использование нестыкующихся на границах подобластей сеток позволяет сформулировать адаптивную технологию, позволяющую строить сетки независимо по подобластям, т.е. параллельно. При этом основная тяжесть параллелизации переносится на создание эффективных параллельных технологий решения возникающих сеточных систем, чему посвящены первые два раздела главы. Формулировка и тестирование адаптивной технологии с использованием нестыкующихся сеток вынссены в третий раздел. Случай декомпозиции области со стыкующимися сетками рассмотрен в разделе 3.4, где предложен новый параллельный метод решения сеточных систем на основе мозаично-скелетонных аппроксимаций плотных матриц. Сужение класса симплициальных адаптивных сеток до иерархических локально сгущающихся сеток позволяет сформулировать в разделе 3.5 еще одну параллельную адаптивную технологию, основанную на многосеточном алгоритме. Используемые на практике консервативные смешаные конечно-элементные аппроксимации порождают системы уравнений, эквивалентные неконформным аппроксимациям, на которые и опирается разработанная технология.

В разделе 3.1 дается макро-гибридная формулировка краевых (диффузионных) задач па нестыкующихся сетках [15], а также общая технология параллельного итерационного решения соответствующих линейных систем. Аппроксимации на нестыкующихся сетках являются неконформными, поскольку используемые конечно-элементные пространства порождают разрывные на интерфейсах решения, которые не принадлежат пространству обобщенных решений (Л). Разобъем область П на тп непересекающихся подобластей Qi регулярной формы и рассмотрим в П краевую задачу

ди

- V ■ рЧи + ей = / в Г2, —- = 0 на дП,

оп

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

Макро-гибридная конечно-элементная аппроксимация этого уравнения приводит к

системе линейных уравнений в седловой форме:

Л,

В,

' А ВТ ' и ' / "

В 0 А 0

0 BJ " щ

Am Bl um

Bm 0 А 0

ся блочно-диагональным переобуславливателем R ■

, для которого необхо-

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

Для построения быстросходящегося итерационного алгоритма можно воспользовать-

" ra О О Лл _

димо определить переобуславливатели в подобластях Йд и на интерфейсе Rx для матриц А и ВА~1ВТ, соответственно: Ra ~ А, ~ ВА~1ВТ.

Переобуславливатель в подобластях Ra выберем блочно-диагональным с блоками

Ri, переобуславливающими матрицы порядка щ. Отметим, что Ai = piAi + где Aj и — матрицы жесткости и масс в подобласти Пусть А'; — переобуслав-

о о

ливатель для матрицы Ai + ¡¿тМ;: c^Ki < Ai + jjMi < с2А'(, где С] > 0, dt — диаметр подобласти регулярной формы Определим матрицу Р; = wi^wf^ через вектор = Q¡е^ (MiWij, wlti) = 1, е; = [1 ... 1]т € R"', а{ = Лемма 3.5. Пусть П^1

е~1Р,, тогда R., ^ А, с константами, зависящими от с\, c-i и не зависящими от pi,£i,di,ni.

Таким образом, задача построения переобуславливателя Ri для оператора —V'pV+г/ в подобласти диаметра di свелась к построению переобуславливателя для оператора —Д-1-1 на растянутой до диаметра 1 сетке f2j, что осуществимо многими способами (см. раздел 1.1). Технология построения переобуславливателя в подобластях Ra обоснована.

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

Пусть матрица ВЙВТ спектрально эквивалентна матрице ВА~1ВТ, однако эффективное решение системы с матрицей ВНВТ невозможно, в отличие от эффективного

• произвольная конформная сетка, Ef < , Rf1 = р.^А',-1 +

умножения на эту матрицу. Определим интерфейсный переобуславливатель Н\ следующим образом (Ю.А.Кузнецов)'.

/?а = внвт (Од - Ц(/л - ъп.-х1внвт)^ ,

что соответствует итерациям ричардсоновского процесса с нулевым начальным приближением. Если в качестве выбрать чебышевские параметры, то число итераций необходимых для Я\ ~ ВЙВТ, должно быть порядка где и — число обусловленности Й^1 В Н Вт. Выбор матриц В и Н базируется на предположении, что наряду с сетками Л* имеются разреженные сетки (необязательно конформные), имеющие тот же след на Г* = д^, что и П , но содержащие гораздо меньше узлов. Построение таких сеток, как правило, не является сложной задачей. Матрица В является аналогом мат рицы В на сетках П^, а матрица Н есть либо факторизованный аналог матрицы А н сетках , либо его переобуславливатель оптимального порядка вычислительной слож ности. В диссертации привены оценки арифметической сложности внутренних итераци: для некоторых наиболее распространненых случаев сеток. Аналогично выводу переобу-славливателя в подобластях воспользуемся следующим утверждением.

Лемма 3.6. Пусть разреженная сетка такова, что для нее выполняется сеточный аналог теоремы о следах, а матрица = К"[ > 0 спектрально эквивалентна

матрице + < Цр, Т\ := [/г, 0] 6 Нп<х"г<. Тогда для дополнения по Шуру

матрицы Л; е П"1*"', 5Г( = Лг< — Лг,/, Лд1 ЛдГ, 6 КПГ|Х"Г', верно

Рх

где РГ( = ш^г-.и^г., и)1,г, = Дег,, г,, ™1.гЛ = 1. еГ( = [1 ... 1]г <= Я"г', А =

| > — интерфейсная матрица масс. Константы спектральной эквивалентности не зависят от р\, £;, й\, щ.

С помощью леммы легко показать спектральную эквивалентность матриц ВА_1ВТ и ВИВ* = £ дГ( + ХГ.А'Г1^) В1-

Легко обратимая матрица Яд = X) Вг< (-гт Рг, + ^'Ц7,1*) ^г введена во внутренние

<=1 V л / •

итерации для устранения зависимости величины и от таких параметров задачи, как

скачки коэффициента диффузии pi, малость коэффициента реакции , большое число

подобластей ш, малый диаметр подобластей

Теорема 3.7. Пусть выполнены условия Леммы 3.6. Тогда число обусловленности и

матрицы {1^ВНВТ не зависит от рх, £i, (1,, и пропорционально = , где

А^г, — собственные значения задачи к,глМг,и>гш, — дополнение по Шуру

сеточной аппроксимации оператора —Д на сетке .

Для случая квази-равномерных сеток с шагом И имеем А„г>Г;/А2,г, ~ <и/к, и для решения системы с матрицей ВНВТ с заданной точностью 0(1) необходимо выполнить Ьх ~ \fdijh внутренних чебышевских итераций, причем Ь\ ~ \J~d~fh не зависит от параметров задачи р<, £;, т, .

Параллелизация внутреннего итерационного процесса аналогична параллелизации решения систем с переобуславливателем Па '■ для этого достаточно распределить разреженные сетки в подобластях по тем же процессорам, что и исходные сетки. Дополнительная сложность привносится необходимостью решать системы с матрицей Яд, что реализуется за счет решения некоторой каркасной задачи собиранием данных (одно число на (подобласть) со всех процессоров на один, решением каркасной системы и распределением данных (одно число на подобласть) с выделенного процессора всем остальным.

Альтернативой внутреннему итерационному процессу может служить Дирихле-Дирихле переобуславливатель. Из Леммы 3.6, примененной к исходной сетке П*, следует, что для интерфейсного блока верно:

т т т ,

ВА'1ВТ = ^В<А~1ВТ = ~ -±_ВГ(РГ(В;С + С,

(=1 ¿=1 '

где <7 = 2 Вт, МГ( - Аи1,А^А1{, ЛГ( - ^пл^/.'^лг, — дополнение по Шуру

о 1

для матрицы А{ = A¡ +

Теорема 3.8. Пусть сетка О;1 такова, что для нее выполняется сеточный аналог теоремы о следах, а матрица

> 0 такова, что спектпр матрицы ЗС пуптадле-жит отрезку [с^сз], 0 < < с^ и пусть

1=1

Тогда ~ ВА~1ВТ. Константы спектральной эквивалентности не зависят от /л, £(, ггр(, т, но зависят от С1,С2.

Вывод Дирихле-Дирихле переобуславливателя 7 для С является эмпирическим. В случае рс = 1, » = 1,..., т:

•Л = ]Г Р^Ви^и (^г, - АглА^Аьъ) "г •=1

^К = ЫосЫЬДОг«,*}. = Вг^шг, В?^ + В^..

Здесь элементы диагональных матриц о>Г; обратны к числу подобластей, которым принадлежит соответствующий узел на интерфейсе, обозначает область, соседнюю с Г2(, с общими гранями ] и , а .Вг^ — вклад грани ] в интерфейсную матрицу £?г( • Факторизация матрицы Fгi недорога, поскольку Рг, — разреженная матрица. Блочная конструкция переобуславливателя обеспечивает его легкую параллелизацию. Наличие матрицы не усложняет параллелизации, поскольку она блочно-диагональна и оперирует исключительно с распределенными по процессорам данными. Рассмотрены также эмпирические обобщения метода на случай разрывных коэффициентов диффузии Pi.

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

Сочетание параллельной итерационной технологии с адаптивной генерацией расчетных сеток независимо по подобластям порождает новый класс параллельных адаптивны? алгоритмов, описанный в разделе 3.3. В дополнение к уже изложенному инструментарию (процедура решения и макро-гибридная дискретизация) предположим, что нам доступен адаптивный генератор симплициальных сеток в расчетной области, а также в подобластях-симплексах.

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

Шаг инициализации. Построить начальную симплициальную сетку, задающую разбиение на т > Р подобластей Распределить (равномерно) подобласти по Р процессорам. В каждой подобласти построить квази-равномерную сетку с примерно одинаковым числом элементов.

Итерационный Шаг.

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

2. Инициализировать переобуславливатель Я, формируя локальные переобуслав-ливатели и оценивая арифметическую нагрузку на каждой подобласти.

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

4. Перераспределить подзадачи по процессорам согласно новому распределению.

5. Решить систему уравнений.

6. Адаптивно перестроить сетку в каждой подобласти независимо друг от друга и перейти к 1.

— собственные значения задачи Зг.^г, = А^М^гиг,. ¿г, — дополнение по Шуру сеточной аппроксимации оператора —А на сетке П^1.

Для случая квази-равномерных сеток с шагом К имеем А„г ,г,/-^2,Г1 ~ и для

решения системы с матрицей ВНВТ с заданной точностью 0(1) необходимо выполнить Ьх ~ у/сЦ/К внутренних чебышевских итераций, причем Ьх ~ ^/ТЦ/Ть не зависит от параметров задачи , т, .

Параллелизация внутреннего итерационного процесса аналогична параллелизацин рс-шения систем с переобуславливателем Р.^: для этого достаточно распределить разреженные сетки в подобластях по тем же процессорам, что и исходные сетки. Дополнительная сложность привносится необходимостью решать системы с матрицей Их, что реализу-гся за счет решения некоторой каркасной задачи собиранием данных (одно число на одобласть) со всех процессоров на один, решением каркасной системы и распределением данных (одно число на подобласть) с выделенного процессора всем остальным.

Альтернативой внутреннему итерационному процессу может служить Дирихле-Дирихле переобуславливатель. Из Леммы 3.6, примененной к исходной сетке следует, что для интерфейсного блока верно:

ВА~*ВТ = ¿дд-'вг = £ ВгА-Х ~ £ + б,

1=1 ¡=1 1=1 Е<а<

„ т —1 я л л л

где & = Вг. (Лг( — ) Вр., Лг, — -^г./.^^М/.г, — дополнение по Шуру

для матрицы A¡ = Ai + ¿гМ¡.

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

> 0 такова, что спектр матрицы 3<3 принадлежит отрезку [сх, Сг], 0 < Сх < сг и пусть

т 1

Тогда Я\ ~ ВА~*ВТ. Константы спектральной эквивалентности не зависят от £(, , пр,, т, но зависят от с^ с2.

Вывод Дирихле-Дирихле переобуславливателя 3 для 6 является эмпирическим. В случае р( = 1, г = 1,..., т:

•А = (^Г, - шг.В^г"1,

>=1

Здесь элементы диагональных матриц обратны к числу подобластей, которым принадлежит соответствующий узел на интерфейсе, Пи обозначает область, соседнюю с Г2;, с общими гранями з и а Вг,^ — вклад грани з в интерфейсную матрицу £?г,. Факторизация матрицы Р^ недорога, поскольку Р^ — разреженная матрица. Блочная конструкция переобуславливателя обеспечивает его легкую параллелизацию. Наличие матрицы РГ( не усложняет параллелизации, поскольку она блочно-диагональна и оперирует исключительно с распределенными по процессорам данными. Рассмотрены также эмпирические обобщения метода на случай разрывных коэффициентов диффузии р,.

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

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

Шаг инициализации. Построить начальную симплициальную сетку, задающую разбиение на тп ~3> Р подобластей П,. Распределить (равномерно) подобласти по Р процессорам. В каждой подобласти построить квази-равномерную сетку с примерно одинаковым числом элементов.

Итерационный Шаг.

4. Перераспределить подзадачи по процессорам согласно новому распределению.

5. Решить систему уравнений.

6. Адаптивно перестроить сетку в каждой подобласти независимо друг от друга и перейти к 1.

Оптимальное распределение подзадач по процессорам осуществляется решением следующей оптимизационной задачи: разбить массив нагрузок Q = {?i,...,9m} па Р подмножеств Pt таким образом, чтобы индекс балансировки min Г j/ max У) ij; был

t=i,...,pi€n .....

как можно ближе к 1. Эта задача может быть эффективно решена в предположении, что m 3> Р. Предложенная технология была успешно протестирована на ряде краевых задач.

В последних разделах третьей главы рассматриваются неконформные аппроксимации на конформных (стыкующихся) сетках. В разделе 3.4 используются макро-гибридные формулировки с неконформными множителями Лагранжа, обеспечивающие поузловую непрерывность сеточного решения и его потока на интерфейсе. С учетом конформности сетки можно показать, что полученное решение удовлетворяет стандартной юнечно-элементной системе па стыкующихся сетках. Для построения интерфейсного пе-реобуславливателя Rx ~ ВА~1ВТ в некоторых случаях удается использовать нормировку сеточных функций в пространстве следов (C.B.Непомнящих) на интерфейсах. Эффективность применения R\ достигается технологией мозаично-скелетонной аппроксимации (В.Е.Тыртышников). Основная идея технологии состоит в аппроксимации плотной матрицы, возникающей из нормализации следа, некоторой разреженной матрицей, которую легко умножать на вектор. Реализация переобуславливателя, использующая в качестве входных данных коэффициенты задачи и интерфейсную сетку, была протестирована на ряде трехмерных краевых задач и показала свою эффективность.

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

1. на входе дается линейная смешаная конечно-элементная система в дисассемблированном (по ячейкам сетки) виде, а также координаты вершин симплексов и граф связности вершин;

2. с использованием геометрических данных система преобразуется локально по ячейкам в гибридную систему за счет добавления множителей Лагранжа на грани симплексов;

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

4. итерационно решается разреженная конденсированная система для множителей Лагранжа;

5. восстанавливаются локально по ячейкам скалярные неизвестные внутри симплексов и потоки на гранях симплексов.

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

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

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

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

В разделе 4.1 формулируется и анализируется параллельный двухуровневый метод Шварца в применении к сингулярно возмущенному уравнению конвекции-диффузии. Предлагаемый алгоритм базируется на том, что возмущения решения этого уравнения, вызванные локальным возмущением краевого условия, распространяются анизотропно. В случае постоянного вектора переноса точечное возмущение распространяется только вниз по потоку, быстро убывая в поперечном направлении и вверх по потоку. На первом уровне метода область разбивается на налегающие полосы (слои) Clk, к = 1,... ,т, перпендикулярные доминирующему направлению переноса, и определяются итерации Шварца в направлении вниз по потоку. На втором уровне каждая подобласть-полоса разевается на налегающие подобласти -квадраты s — 1,...,р, и задаются несколько „Итераций Шварца между квадратами с четными и нечетными индексами я. Естественный параллелизм методу дает именно второй уровень, поскольку каждая подобласть-квадрат может обрабатываться своим процессором. Если вектор переноса сильно отклоняется от заданного направления (например, вследствие завихрений), разбиение первого уровня должно быть адаптировано к поведению вектора. Доказана универсальность метода по отношению к большим числам Пекле.

Рассмотрим в области U = (0;1) х (0; 1) двумерное уравнение конвекции-диффузии с оператором С — —еД + д—■ (е = const 6 (0; 1)) с естественными и главными краевыми условиями на Tjv и ЗП\Глг, соответственно, и аппроксимирующую систему разностных уравнений на квадратной сетке Clhl ft = l/n»£:

Uij = fij, Xij £ U U VN, иц — g^, x{j edU\ TN, xi} = (ih,jh), i,j = 0,..., n.

Сформулируем двухуровневый метод Шварца для случая = ("t.i'i) х (0,1), разбитых последовательностью целых чисел 0 = n'1' < n^1' = nj2' < ... < т4''' = п на ilk, —

(ак-,Ьк) X (n^hih^h), s четное, .(>) w (я)

; , ! .), и,\ . где п\' = max{0;n\' - dh}, «V = +

(o.k',bk) x (nj h-,n2 Л), s нечетное

dh}, dh = min (6* — at+i) — минимальная ширина налегания полос Q/,. Пусть задано

l<fc<m—1

текущее приближение u'j1 к сеточному решению и¡¡. Новое приближение u'j находится с помощью двухуровневого метода Шварца:

Алгоритм. Пусть заданы целые Тк > 0- Для к = 1,...,т найти сужение сеточной функции u\j на (ак\ aj.+i] х (0; 1) следующими внутренними итерациями:

1. Положить 1 — 0. Для всех четных s решить подзадачи -4^ „iiy' = fij.

2. Повторять пока ( I <Тк )

Положить I = I + 1.

{

Если I нечетно, то

для всех нечетных s решить подзадачи ЛЛ'^й^ — fij иначе

для всех четных s решить подзадачи .Mj^-ûy' = /у Конец цикла

3. Определить сеточную функцию и{• на [at, at+i] X (0; 1) через u\j = ûy\ х € К, iifc+i] X [n'/'fc; т4"'л], s = 1,..., р.

Здесь через Ai^û^f = /у обозначено сужение на Як, исходной разностной задачи с главным краевым условием и.у там, где «у уже известно (xi = ак), или ujj1 в противном случае. Отличие подзадач = /у от — fij заключается в том, что там

где известно итерационное приближение ûy*1' к Uy со стороны соседних подобласте! ii краевые условия с заменяются на ûy±l'.

Для случая s <S. h сходимость метода с Tt = 1, к — 1,..., т, устанавливается следствием более общего утверждения (Теорема 4.1)

Следствие 4.5. Пусть s <S h, Тк = 1, к = l,...,m, /у > 0, jy > 0, и задана неотрицательная сеточная функция и?"1 на Ù такая, что тахщ, |"у — —

Тогда

max |uy — Uy| < Zmq~[d,'/h max |uy — ujj1!, qx = -.

Таким образом, для сингулярно возмущенного оператора конвекции-диффузии с постоянным конвективным полем геометрическая скорость сходимости метода Шварца обеспечивается минимальными перекрытиями между подобластями и всего лишь одной внутренней итерацией в полосах Пк■

В случае переменного конвективного поля с доминирующим направлением геометрическая скорость сходимости метода Шварца обеспечивается за счет адаптивного разбиения первого уровня на подобласти, точного решения подзадач в подобластях-полосах с противотоком, минимальными перекрытиями между подобластями-полосами и всего лишь одной внутренней итерацией в полосах iijt без противотока, где налегание подобластей второго уровня Çlk,i должно быть порядка ширины iIk .

Оператор задачи с переменным конвективным полем имеет вид С = —E'Д-|-.B1(2:l,£2)-g^ + В2(:Г],^2) - с дополнительным предположением, что функции В1 и В2 достаточно гладкие, не превышающие по модулю 1, и B'Çl.Xj) > 0. Разбиение Q на полосы Пк должно быть адаптировано (что может быть реализовано автоматически) таким образом, что вектор переноса В = (ВВ2)Т в О^ удовлетворяет одному из двух условий: либо B}j > В о > 0, -Bfj < R ■ В0 в П it, с некоторой константой 0 < R < 1, либо

B}j > Во > 0, в fit \ (at;at+i) х (0; 1). В первом случае индекс принадлежит множеству V,, во втором — множеству "Pv. Формулировка двухуровневого метода Шварца ничем не отличается от случая постоянного поля переноса. В случае переменного поля мы определим последовательность Т,t следующим образом:

i Т, к I \ 1, к <

Тк = { " '' € ^

где Т — некоторое большое целое число. Предположение Т 3> 1, означает, что задачи в подобластях $1к, к £ Р,,, решаются точно. Упрощенное для случая £ < Л утверждение " сходимости метода выглядит следующим образом:

Теорема 4-6. Пусть Г> 1, /у > 0, > 0, разбиение па С1к адаптировано

В, и задана неотрицательная сеточная функция на Г2 такая, что тахлЛ

и-~Ч < 1. Если

, . ,, 1пСш . bk-ak+1 1п2#Р„ Во/,

<гЛ = min (б* - dt+i) > h ■ --, min--- > H--;-, <?:i = H--,

i<*<m-iv ln<j3 (сея, /i ln<j3 £

ter.

[n^ - г4")Л = К" - n\'yi > max (bt - а*) +

ie p.

dhlnqa h ln

1 + ^f* 94 = -%-r,

1 + Я- ^

¿.C

п* — —

с константой <3 < 1, не зависящей от £, к.

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

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

В разделе 4.3 представлено несколько адаптивных итерационных технологий, ускоряющих время решения последовательности систем, порождаемых в процессе дискретизации по времени. Эти технологии легко параллелизуются и не используют свойства расчетной сетки: их основной целью является построение хорошего начального приближения для текущей системы на основе информации, полученной при решении предыдущих систем. Рассматриваются три типа последовательностей систем уравнений: (а) линейные системы с одной квадратной матрицей А и разными правыми частями Ьк\ (б) линейные системы с разными, но близкими квадратными матрицами Ак и разными правыми частями Ьк\ (в) нелинейные системы.

В первом случае используется обобщенный метод сопряженных невязок, где помимо ЛгЛ-ортогональных крыловских векторов К. = {Pj}JLn известны векторы {Apj}^Li. Поскольку проекция Ьк на пространство АК есть j>2 (Ар^Ар)^'' то проекция решен^

хк — А~1Ьк на К. есть хк = (ар'^Ар1)Рз- Таким образом, проекция неизвестного решени

на накопленное крыловское подпространство может быть легко вычислена посредство... к скалярных произведений, а накопление подпространств К и АК может быть продолжено для разных векторов Ьк, если начальное приближение для каждого последующего решения есть хк. Эта стратегия приводит к 1.5-2-кратному ускорению расчета больших жестких линейных систем, порождаемых проекционным алгоритмом для решения нестационарных уравнений Навье-Стокса.

Во втором случае используется обобщенный метод минимальных невязок GMRES и накапливаются базисы крыловских пространств Vmt и данные о проекции Ак на Vmi, матрицы Gmt и IImk,mt. Начальное приближение хк для fc-й системы вычисляется последовательностью проектирований = if, — — 6'), i = к —I,...,к —I, где Xq~1 = 0. Число накапливаемых данных I ограничено не только емкостью компьютерной памяти, но и расхождением между собственными векторами и значениями матриц Ак~1 и А1"1. Эта стратегия приводит к 15-процентному ускорению расчета больших жестких линейных систем, порождаемых проекционным алгоритмом для решения уравнения тепло-массопереноса в простейшем химическом реакторе.

В третьем случае рассматривается применение метода Ньютона для полностью неявных дискретизаций нелинейных нестационарных краевых задач. Физически мотивированный выбор начального приближения — использование решения с предыдущего временного шага. Для нахождения лучшего приближения воспользуемся собственным ортогональным разложением (POD). POD порождает оптимальные аппроксимации малых размерностей для наборов данных, генерируя ортонормированный базис для представления данных с помощью метода наименьших квадратов. Используя галеркинскую про-

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

В разделе 4.4 излагаются параллельные технологии эффективного решения систем уравнений многофазной фильтрации как многоуровневого, так и блочного типов. В случае неявных дискретизаций многофазных многокомпонентных течений, порождаемые методом Ньютона линейные системы содержат уравнения, соответствующие различным процессам, происходящим в разных масштабах, поэтому прямой перенос современных технологий (многосеточные методы, методы декомпозиции), разработанных для эффективного решения эллиптических уравнений, невозможен. Рассмотренный метод алгебраического выделения уравнения для давления имеет два потенциальных преимущества: во-первых, можно строить хороший переобуславливатель только для самых жестких блоков, переобуславливая другие блоки вычислительно дешевыми технологиями; во-вторых, этот подход не зависит от структуры сетки, поскольку переобуславливатели для блоков могут опираться только на их алгебраическую структуру. Второй предлагаемый метод, многосеточный метод с ускоренным огрублением SCMG, использует иерархичность трехмерной прямоугольной сетки и вертикальное агрегирование неизвестных. Несмотря на структурное подобие SCMG и стандартного V-цикла геометрического многосеточного метода, они опираются на разные итерационные механизмы. SCMG оказался эффективным летодом решения задач с сильно меняющимися коэффициентами, несмотря на простей-пие операторы интерполирования данных между уровнями сетки.

Основные результаты четвертой главы состоят в следующем: (а) параллельный двухуровневый метод Шварца в применении к сингулярно возмущенному уравнению конвекции-диффузии обобщается и обосновывается на случай произвольных полей вектора скорости с доминирующим направлением конвекции; (б) рассмотрены параллельные трехмерные приложения разработанной технологии; (в) рассмотрены и протестированы три параллельные стратегии адаптивного выбора начального приближения при итераци-

онном решении последовательностей линейных систем с одной матрицей, разными матрицами, и нелинейных систем; (г) разработаны и исследованы две параллельные технологии решения линейных систем, возникающих при полностью неявных аппроксимациях уравнений многофазной фильтрации.

Основные результаты диссертации

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

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

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

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

4. рассмотрены теоретические и практические вопросы управления адаптацией в рамках вышеупомянутой адаптивной технологии; доказаны оценки сверху для интерполяционных ошибок на сетках с управляемой адаптацией;

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

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

рующим направлением поля скорости; доказана универсальность метода по отношению к большим числам Пекле; рассмотрен ряд трехмерных приложений метода;

7. предложены, реализованы и протестированы параллельные технологии адаптивного выбора начального приближения при итерационном решении последовательностей линейных систем с одной матрицей, разными матрицами, и нелинейных систем;

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

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

1. Василевский Ю.В. Методы решения краевых задач с использованием нсстыкую-щихся сеток. // Тр. Матем. центра им. Н.И.Лобачевского, Т.2. Итерационные методы решения линейных и нелинейных сеточных задач. - Казань: УНИПРЕСС, 1999, С.94-121.

2. Василевский Ю.В., Агузал А. Объединенный ассимптотический анализ интерполяционных ошибок на оптимальных сетках. // ДАН, 2005, Т.405, No.3, С.1-4.

3. Василевский Ю.В., Липников К.Н. Адаптивный алгоритм построения квазиоптимальных сеток.// ЖВМ и МФ, 1999, Т.39, No.9, С.1532-1551.

4. Василевский Ю.В., Липников К.Н. Оптимальные триангуляции: существование, аппроксимация и двойное дифференцирование Pi конечно-элементных функций. // ЖВМ и МФ, 2003, Т.43, No.6, C.8CG-874.

5. Василевский Ю.В., Липников К.Н. Оценки ошибки для управляемых адаптивных алгоритмов на основе восстановления гессиана. // ЖВМ и МФ, 2005, Т.45, No. 8, С.1424-1434.

6. Василевский Ю.В., Капранов С.А. Параллельное моделирование особенностей кровотока в окрестности кава-фильтра с захваченным тромбом. // Математическое моделирование, 2005, Т.17, No.ll, С.3-15.

7. Agouzal A., Lipnikov К., Vassilevski Yu. Adaptive generation of quasi-optimal tetrahedral meshes. // East-West J. Numer. Math., 1999, V.7, P.223-244.

8. Agouzal A., Vassilevski Yu. On a discrete Hessian recovery for Pi finite elements.// J. Numer. Math., 2002, V.10, P.l-12.

9. Chugunov V., Vassilevski Yu. Parallel multilevel data structures for a non-conforming finite element problem on unstructured meshes.// Russ.J.Numer.Anal.Math.Modelling, 2003, V.18, No.l, P.l-11.

10. Chugunov V., Svyatski D., Tyrtyshnikov E., Vassilevski Yu. Parallel iterative multilevel solution of mixed finite element systems for scalar equations.// Concurrency and Computation: Practice and Experience, 2006, V.18, No.5, P.501-518.

11. Dyadechko V., Iliash Yu., Vassilevski Yu. Structuring preconditioners for unstructured meshes. // Russ.J. Numer.Anal. Math.Modelling, 1996, V.ll, No.2, P.139-154.

12. Dyadechko V., Lipnikov K., Vassilevski Yu. Hessian based anisotropic mesh adaptation in domains with discrete boundaries. // Russ. J.Numer. Anal.Math.Modelling, 2005, V.20, No.4, P.391-402.

13. Garbey M., Kuznetsov Yu., Vassilevski Yu. Parallel Schwarz method for a convection-diffusion problem.// SIAM J.Sci.Comp., 2000, V.22, No.3, P.891-916.

14. Garbey M., Vassilevski Yu. A parallel solver for unsteady incompressible 3D Navier-Stokes equations.// Parallel Computing, 2001, V.27, No.4, P.363-389.

15. Hoppe R., Iliash Yu., Kuznetsov Yu., Vassilevski Yu., Wohlmuth B. Analysis and parallel implementation of adaptive mortar element methods. // East-West J. Numer. Math., 1998, V.6, No.3, P.223-248.

16. Lacroix S., Vassilevski Yu., Wheeler M., Wheeler J., Iterative solution methods f" modeling multiphase flow in porous media fully implicitly.// SIAM J.Sci.Comp., 2003, V.25, No.3, P.905-926.

.17. Lipnikov K., Vassilevski Yu. Parallel adaptive solution of 3D boundary value problems by Hessian recovery.// Comp.Methods Appl.Mech.Engnr., 2003, V.192, No.11-12, P.1495-1513.

18. Lipnikov K., Vassilevski Yu. Parallel adaptive solution of the Stokes and Oseen problems oil unstructured 3D meshes. // Proceedings of Int.Conf. Parallel CFD 2003, Elsevier, 2004, P. 153-162.

19. Lipnikov K., Vassilevski Yu. On control of adaptation in parallel inesh generation. // Engrg. with Comput., 2004, V.20, No.3, P.193-201.

20. Tromeur-Dervout D., Vassilevski Yu. POD acceleration of fully implicit solvers for unsteady nonlinear problems and its GRID applications. // Advances in Engineering Software, 2006, to appear.

21. Tyrtyshnikov E., Vassilevski Yu. A mosaic preconditioner for a dual schur complement.// Numerical Mathematics and Advanced Applications, Proc. of ENUMATH'01. -Milano: SpringerVerlag Italia, 2003, P.867-880.

22. Vassilevski Yu. A hybrid domain decomposition method based on aggregation.// Numer. Linear Algebra Appl., 2004, V.ll, P.327-341.

23. Vassilevski Yu. A parallel interface preconditioner for the mortar element method in case of jumping coefficients.// Domain Decomposition Methods in Sciences and Engineering, DDM.org, 2001, P.231-240.

24. Vassilevski Yu. A parallel CG solver based on domain decomposition and non-smooth aggregation. // Conjugate gradient algorithms and finite element methods (Proceedings of Int.Conf. 50 years of CG), -Berlin-Heidelberg: Springer-Verlag, 2004, P.93-102.

Изд.лиц. ИД N 03991 от 12.02.2001. Компьютерный набор Подписано в печать 27.02.2006. Усл. пет. л. 2,0. Тираж 80 экз. Институт вычислительной математики РАН 119991 ГСП-1, г.Москва, ул.Губкина 8.

Содержание диссертации автор исследовательской работы: доктора физико-математических наук, Василевский, Юрий Викторович

Введение

Глава 1 Блочные параллельные методы для решения неструктурированных систем

1.1 Последовательные технологии решения неструктурированных систем

1.1.1 Методы факторизации

1.1.2 Методы неполной факторизации.

1.1.3 Алгебраические многосеточные методы.

1.1.4 Структурированные переобуславливатели для неструктурированных сеток

1.2 Решение неструктурированных систем с матрицами тензорного ранга

1.3 Блочные переобуславливатели и переупорядочивание для задач с анизотропными коэффициентами

1.4 Параллельный метод агрегирования и его применения.

1.4.1 Метод агрегирования.

1.4.2 Параллелизация метода агрегирования

1.4.3 Диффузиоипые задачи с гетерогенными коэффициентами и/или большим числом подобластей

1.4.4 Параллелизация на адаптивных неструктурированных сетках . 65 Выводы главы 1.

Глава 2 Параллельные адаптивные технологии на симплициальных анизотропных сетках

2.1 Оптимальные и квази-оитимальные сетки и их свойства.

2.1.1 Оптимальные сетки.

2.1.2 Квази-оптимальные сетки

2.2 Последовательный адаптивный алгоритм построения квази-оптимальных сеток.

2.2.1 Восполнение сеточного гессиана

2.2.2 Генерация сеток, квази-равномерных в заданной метрике

2.2.3 Адаптивный алгоритм.

2.3 Управление адаптацией.

2.4 Особенности адаптации для областей с дискретной границей.

2.5 Параллельный адаптивпый алгоритм для приближенного решения краевых задач .ИЗ

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

2.5.2 Параллельное адаптивное решение трехмерных краевых задач

2.5.3 Влияние управления адаптацией

Выводы главы 2.

Глава 3 Блочные параллельные методы решения неконформных конечноэлементных систем

3.1 Макро-гибридная постановка па песты кующихся сетках.

3.2 Интерфейсные переобуславливатели

3.2.1 Переобуславливатели с внутренним итерационным процессом для множителей Лагранжа.

3.2.2 Дирихле-Дирихле переобуславливатель.

3.2.3 Численные эксперименты с параллельными решателями.

3.3 Параллельные адаптивные технологии на нестыкующихся сетках

3.4 Мозаично-скелетонные переобуславливатели для множителей Лагранжа

3.5 Параллельный многосеточный метод для неконформных конечных элементов 165 Выводы главы 3.

Глава 4 Блочные параллельные технологии на трехмерных прямоугольных сетках

4.1 Параллельный двухуровневый метод Шварца для сингулярно возмущенного уравнения конвекциидиффузии

4.1.1 Двухуровневый метод Шварца для постоянного вектора переноса

4.1.2 Двухуровневый метод Шварца для переменного вектора переноса

4.1.3 Адаптивный алгоритм разбиения для случая переменного вектора переноса.

4.2 Некоторые приложения параллельного двухуровневого метода Шварца . . 193 4.2.1 Трехмерное уравнение конвекции-диффузии.

4.2.2 Трехмерное уравнение конвекции-диффузии-реакции

4.2.3 Проекционный алгоритм для нестационарных уравнений Навье-Стокса

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

4.3 Адаптивные итерационные технологии для последовательности систем

4.3.1 Последовательность систем с одной матрицей и разными правыми частями.

4.3.2 Последовательность систем с разными матрицами и правыми частями

4.3.3 Последовательность нелинейных систем.

4.4 Блочные параллельные методы для систем уравнений многофазной фильтрации

4.4.1 Уравнения многофазной фильтрации, их дискретизация и линеаризация

4.4.2 Многосеточный метод с ускоренным разгрублением

4.4.3 Алгебраическое выделение уравнения для давления.

Выводы главы 4.

Введение диссертация по математике, на тему "Параллельные технологии решения краевых задач"

Настоящая диссертация посвящена рассмотрению некоторых современных вычислительных технологий нахождения приближенного решения краевых задач с эллиптическими операторами второго порядка. Под вычислительными технологиями здесь понимается совокупность алгоритмов, структур данных, расчетных методик и программных реализаций для решения вычислительных задач на вычислительных системах. Главной целью вычислительных технологий является получение ответа в нужной форме, к заданному сроку и при известных ограничениях на доступ к машинным и человеческим ресурсам [8] . Основные черты рассматриваемых здесь технологий — параллельность и эффективность.

Параллельность вычислительных технологий вызвана необходимостью решать задачи настолько большой размерности, что даже новейшие однопроцессорные компьютеры не в состоянии обеспечить требуемый расчет. Помимо практических ограничений на быстродействие процессора и размеры компьютерной памяти, немаловажно учитывать и экономический аспект вычислений: очевидна нелинейная зависимость цены компьютера от объема памяти и реального быстродействия процессора. В силу этого параллелизация вычислений, во-первых, позволяет решать гораздо большие задачи, и во-вторых, обеспечивает примерно линейную зависимость цены параллельного компьютера с распределенной памятью от его быстродействия и объема памяти. Именно последнее обстоятельство объясняет широкую распространенность параллельных компьютеров с распределенной памятью и естественный запрос на параллельные технологии, позволяющие использовать такие компьютеры. Распределенная память подразумевает разбиение данных на блоки, обрабатываемые каждым процессором, поэтому блочность алгоритмов характерна для большинства параллельных методов и технологий [18].

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

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

Общепринятым подходом при разработке современных технологий решений краевых задач математической физики является их замена сеточными уравнениями, обеспечивающая приближенное решение этих задач. Одной из основных составляющих технологии построения систем сеточных уравнений является генерация расчетной сетки. Адаптивно построенные сетки позволяют достигать желаемой точности при приближении исходной краевой задачи на меньшем количестве узлов и элементов сетки. Наиболее широкое освещение вопросов адаптивного приближенного решения краевых задач можно найти в книгах Р.Ферфурта [260], И.Бабушки и Т.Струболиса [59], хотя некоторые ключевые моменты технологии были рассмотрены уже в монографии Л.А.Оганесяна и JI.А.Руховца [33]. Вопросы генерации адаптивных сеток рассматривались в монографиях П.Джорджа и Х.Буручаки [143], В.Д.Лизейкипа [193], Г.Кэри [94], а также в сборнике [274]. Системы сеточных уравнений могут решаться как прямыми, так и итерационными методами, однако применение прямых методов ограничено лишь линейными системами с матрицами умеренных порядков или специального вида. Хорошее представление о состоянии теории и практики итерационных методов решения сеточных уравнений можно получить из книг Д.Янга [272], Г.И.Марчука [32], Г.И.Марчука и В.И.Лебедева [30], А.А.Самарского и Е.С.Николаева [37], Е.Г.Дьяконова [22], В.Хакбуша [157], И.Саада [234], Дж.Деммеля [20].

Перейдем к общему обоснованию предлагаемых в диссертационной работе технологических решений.

Использование адаптивных сеток для построения сеточных уравнений зачастую подразумевает решение неструктурированных систем, графы матриц которых не обладают никакой структурой. Аппроксимации па структурированных сетках позволяют использовать эффективные методы решения линейных систем. Так, дискретизации на прямоугольных сетках предполагают лексикографическое упорядочивание узлов и связанных с ними неизвестных в уравнении, что позволяет использовать быстрые прямые методы для решения соответствующих систем [93, 62, 259, 28, 180, 231]. Если сетка представляет собой подмножество ячеек прямоугольной сетки (как в методе фиктивных компонент), то возможно такое упорядочивание узлов, что быстрые прямые методы применимы в качестве переобуславливателей [197, 6]. Использование иерархических локально сгущающихся сеток позволяет использовать многосеточные технологии [40, 34, 158], основанные на интерполяции данных между сетками разных уровней. Отметим усиленные многосеточные методы [159, 183], обеспечивающие высокую скорость сходимости на частично неструктурированных сетках, получаемых двух-, трех-уровневым равномерным измельчением достатачно мелкой неструктурированной сетки. Во всех упомянутых случаях, по крайней мере для модельных задач, удается построить и обосновать либо прямой, либо быстро сходящийся итерационный метод, в котором арифметическая работа одной итерации (или всего метода) пропорциональна числу неизвестных, умноженному на полилогарифмический множитель (в случае прямых методов), что является хорошим асимптотическим результатом.

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

Неструктурированные сетки как симплициальные разбиения областей привлекательны в инженерных приложениях по разным причинам. Сложность расчетной области является одной из ключевых: ни прямоугольная, ни симплициальная грубая сетка, порождающая иерархическую последовательность сеток, зачастую не в состоянии аппроксимировать все характерные особенности границы, обладая приемлемым количеством сеточных узлов. Именно этим объясняется широкое использование генераторов неструктурированных триангуляций в двумерном случае, и тетраэдризаций в трехмерном случае. Второй причиной привлекательности неструктурированных сеток является их гибкость в адаптации. Так, например, возможности адаптации прямоугольных сеток в рамках метода фиктивных областей/компонент исчерпываются сгущением узлов вдоль координатных осей и адаптацией криволинейной границы со вторым порядком точности. Адаптивное построение анизотропных сеток на основе преобразования координат узлов прямоугольной сетки [23, 193, 274] с сохранением связности между узлами также является востребованной технологией построения анизотропных сеток, которая, однако, обладает рядом практических ограничений, особенно в трехмерном случае. Возможности адаптации локально сгущающихся иерархических триангуляций и тетраэдризаций [143], [94], [260] ограничены регулярной формой симплексов: единственная возможность управлять формой элементов — поддерживать регулярность их формы на уровне, заданном начальной неструктурированной сеткой, поэтому подстраивание сеток иод особенности решения имеет смысл для решений с изотропными особенностями. Управление вытянутостью ячеек, или анизотропностью сетки, представляется возможным только для полностью неструктурированных сеток.

Важность анизотропной адаптивности проявляется в краевых задачах с анизотропным решением, например, благодаря наличию погранслоя или входящего двугранного угла [57, 41]. Неструктурированные симплициальные (треугольные и тетраэдральные) конформные (с ячейками, соседствующими по целой общей грани, ребру, или узлу) сетки являются, на наш взгляд, простейшими сетками, способными обеспечить адекватную адаптацию к наиболее широкому классу решений краевых задач. Дальнейшие расширения класса сеток (неконформность, ячейки-полиэдры) не приведут к более тонким оценкам ошибок, хотя и обеспечат большую технологичность аппроксимаций. Адаптивные свойства класса неструктурированных симплициальных конформных сеток рассматриваются во второй главе диссертации. Для этого вводится и исследуется оптимальная сетка как сетка, минимизирующая норму ошибки заданного оператора проектирования среди всех сеток класса с ограниченным сверху количеством элементов. Там же рассмотрены свойства и технологии (в том числе параллельные) построения аппроксимаций оптимальных сеток и некоторые сопутствующие вопросы адаптивной генерации подобных аппроксимаций.

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

Аппроксимации эллиптических уравнений на пестыкующихся сетках порождают линейные системы с седловыми матрицами. Эффективное параллельное итерационное решение подобных систем является самым сложным элементом технологии приближенного адаптивного решения краевых задач на нестыкующихся сетках [9]. В третьей главе разработаны и обоснованы два типа подобных решателей и предложен и протестирован прообраз соответствующей паралельной технологии. С математической точки зрения, ключевыми для использования нестыкующихся сеток являются условия сшивки на границе между подобластями. Использование упрощенных аналогов этих условий для стыкующихся сеток позволяет сформулировать и обосновать новый декомпозиционный алгоритм для конформных сеток и добавить в предлагаемую параллельную технологию дополнительный метод итерационного решения возникающих сеточных систем.

За последнее десятилетие смешанный метод конечных элементов приобрел большое прикладное значение, поскольку он обеспечивает консервативные аппроксимации эллиптических уравнений и одновременное приближение скалярного решения и его потоков [88]. Матрица сеточного оператора является седловой, однако система может быть эффективно и параллельно преобразована статической конденсацией к системе с симметричной (в случае диффузионного оператора) и положительно определенной матрицей, которая для смешанных конечных элементов Равиара-Тома наименьшего порядка [225] аналогична неконформной конечно-элементной аппроксимации Крузэ-Равиара [108]. Параллельная реализация [102] многосеточного метода для конечно-элементных аппроксимаций Крузэ-Равиара на локально сгущающихся иерархических симплициальных сетках позволила построить эффективную параллельную адаптивную технологию [103] решения смешанных конечно-элементных систем. Таким образом, параллельные адаптивные технологии, изложенные в третьей главе, позволяют эффективно применять некоторые практически значимые аппроксимации с использованием нестыкующихся сеток или консервативных смешанных конечных элементов.

Несмотря на гибкость в адаптации и общность аппроксимаций, предоставляемые неструктурированными сетками, многие приложения до сих пор базируются на использовании прямоугольных сеток. Помимо исторических причин, для этого существуют и технологические причины. Действительно, простота реализации, простота аппроксимаций повышенного порядка точности, наличие свойства суперсходимости при аппроксимации гладких решений, минимизация отношения числа сеточных ячеек и числа сеточных узлов, безусловно, могут обеспечить некоторые технологические преимущества прямоугольным сеткам. В четвертой главе рассматривается ряд приложений в решении краевых задач с использованием прямоугольных сеток. Все технологии, предлагаемые в этой главе, не ограничены классом прямоугольных сеток, их изложение в этой главе диктуется либо простотой изложения и анализа на прямоугольных сетках, либо органичностью их представления в рамках определенных приложений. Так, параллельный двухуровневый метод Шварца для сингулярно возмущенного уравнения конвекции-диффузии с переменным конвективным полем проанализирован для конечно-разностных аппроксимаций [137], однако его реализация возможна и для конечно-элементных аппроксимаций на неструктурированных сетках и даже пестыкующихся сетках [162]. Естественным приложением этого метода является параллельное решение нестационарных уравнений Навье-Стокса [140] с помощью проекционного метода [101, 242, 42, 198, 154], в котором необходимо решать сеточный аналог уравнения Пуассона на каждом временном шаге. Это, в свою очередь, приводит к необходимости эффективного решения последовательности больших жестких систем с разными правыми частями. Помимо использования параллельного метода фиктивных компонент, наиболее эффективного на прямоугольных сетках, предлагается использовать адаптивные механизмы выбора начального приближения для итерационных методов. Будучи вычислительно недорогими, технологически простыми и легко па-раллелезуемыми, эти механизмы позволяют существенно сократить число итераций и, следовательно, вычислительную работу. В диссертации рассматриваются три адаптивных метода выбора начального приближения в рамках приложений на прямоугольных сетках, хотя они являются не зависимыми от типа используемых сеток и аппроксимаций.

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

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

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

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

Методы первой группы используют минимальную информацию об элементах матрицы и правой части. К ним относятся прежде всего прямые методы факторизации разреженных матриц, современные реализации которых обсуждаются в [55, 109, 155, 113], итерационные методы неполной факторизации [24],[58],[234], и итерационные алгебраические многосеточные методы [106, 240]. Каждый из методов обладает рядом достоинств и недостатков: методы факторизации фактически не зависят от значений матричных элементов, однако имеют неудовлетворительную асимптотику арифметической работы на этапе инициализации, особенно для матриц, порожденных трехмерными сетками, и большие требования на компьютерную память, налагающие ограничения на порядок матриц; методы неполной факторизации не обеспечивают высокой скорости сходимости, не зависящей от порядка матриц, однако быстро инициализируются; алгебраические многосеточные методы, показывая высокую скорость сходимости и линейную арифметическую сложность как инициализации, так и решения, опираются на определенные свойства исходных матриц, которые не всегда выполняются.

Методы второй группы ("серый ящик") используют дополнительную информацию для построения переобуславливателя. Так, метод фиктивного пространства [205, 206] нуждается в данных о сетке, порождающей сеточный оператор, и соответствующих краевых условиях, при этом порождает спектрально эквивалентные переобуславливатели с временами инициализации и решения, пропорциональными порядку матриц [121].

Можно отметить следующие тенденции, проявляющиеся в последние годы: методы неполной факторизации строят более аккуратное разложение на верхний и нижний треугольные множители [172] за счет более дорогой процедуры инициализации, однако лучшей скорости сходимости; в то же время алгебраические многосеточные методы могут нуждаться в дополнительной информации, например, предоставление матрицы в дисассемблированном виде [87, 96].

Результаты экспериментальных исследований автора [14, 121], проведенных для различных последовательностей матриц увеличивающегося порядка, дали ответы на некоторые практические вопросы по перечисленным выше классам методов. Для нескольких реализаций методов точной факторизации экспериментально выведены асимптотические оценки сложности процедуры факторизации и последующего решения, а также практические ограничения на компьютерную память, для матриц с наиболее популярными структурами разреженности. Это дает возможность оценить целесообразность применения той или другой реализации метода для конкретной задачи. Показано (см. раздел 1.1), что даже наиболее точные процедуры неполной факторизации [172] чувствительны к измельчению и структуре сетки, в то время как метод фиктивного пространства в сочетании с многоуровневой технологией обеспечивает линейную арифмеческую сложность как процедуры инициализации, так и последующего решения.

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

Матрицы неструктурированных систем, порождаемых призматическими сетками, обладают тензорным рангом 2, т.е. нредставимы в виде А\®МъъЛ-А?а ■ Здесь Л23, М23 (Л1, Mi) — матрицы жесткости и масс (лампировапная) на двумерном (соотв.,одномерном) сечении призматической сетки. Для решения таких систем автором предлагается новый метод, обладающий почти линейной арифметической сложностью решения и невысокой сложностью инициализации. Метод представляет собой обобщение быстрого прямого метода [175, 259, 28, 43, 180, 232, 231] на случай неструктурированных матриц Л23 с двумерным графом и использует явный вид матриц Л1; Л23, М\, М23. Смысл обобщения — замена прямого решения двумерных систем с матрицами Л23 + ЛМ23 с произвольным параметром Л на приближенное итерационное решение с нереобуславливателем, базирующемся на факторизации одной из трех матриц Л23+ЮМ23, Л2з+150М2з, Л2з + 1500М2з-При этом результирующий переобуславливающий оператор будет нелинейным, и метод из прямого становится быстро сходящимся итерационным (FGMRES [235]) с нелинейным нереобуславливателем, арифметическая цена применения которого (на практике) пропорциональна порядку матрицы с полилогарифмическим множителем.

Для решения неструктурированных анизотропных систем, порождаемых симплици-альными сетками и анизотропным коэффициентом диффузии, предлагается метод переупорядочивания неизвестных и разбиения матрицы на блоки, с помощью которого простейшие блочные переобуславливатели (блочные методы Якоби, Гаусса-Зейделя) становятся весьма эффективными. Алгоритм, разработанный и реализованный С.А.Горейновым с частичным участием автора, является вариантом жадного алгоритма с ограничениями, перекликающегося с алгоритмами блочного переупорядочивания PABLO [210] и TPABLO [100]. Дополнительные данные, требуемые алгоритмом переуиорядочивания и формирования блоков, — матрица весов, соответствующая ненаправленному графу матрицы системы и задающая анизотропию системы. В случае оператора диффузии в качестве матрицы весов можно взять модули внедиагональных элементов матрицы жесткости. Обладая линейной арифметической сложностью, предлагаемый алгоритм обеспечивает приемлемые скорости сходимости для матриц, порождаемых как сильной, так и слабой анизотропией диффузионного тензора, на сильно сгущающихся неструктурированных двумерных и трехмерных сетках.

Для параллельного решения неструктурированных систем, порождаемых эллиптическим оператором с гетерогенным коэффициентом диффузии на неструктурированных сетках, автором предложен декомпозиционный метод агрегирования [255, 258]. Одной из основных целей декомпозиционных алгоритмов [31, 220, 238, 275] является следующая: при заданном разбиении расчетной области на подобласти и заданных в подобластях решателях или преобуславливателях, решить общую задачу быстросходящимся итерационным методом, легко реализуемым, легко параллелезуемым, и с низкими накладными арифметическими расходами. Желательно, чтобы скорость сходимости итераций не зависела от шага сетки (h), диаметра подобластей (Я), и других параметров задачи, например, скачков коэффициента диффузии (р).

Существуют методы [208], удовлетворяющие иочти всем требованиям, однако они сложны в реализации и имеют определенные ограничения (например, на прохождение границы между подобластями и линиями скачка коэффициента). В силу этого, разумным представляется поиск методов, удовлетворяющих не всем условиям, но желаемым. Например, методы [119,120,195,196, 221] демонстрируют слабую полилогарифмическую зависимость числа обусловленности переобусловленной матрицы жесткости от шага сетки, однако требуют точного решения подзадач, что повышает арифметические затраты, а метод [82] легко параллелезуем, использует в подобластях лишь переобуславливате-ли, однако порождает более сильную зависимость от шага сетки (H/h). Вышеупомянутые методы базируются на разбиении на непересекающиеся подобласти, причем границы раздела подобластей (интерфейсы) предполагаются гладкими и совпадающими с линиями или поверхностями скачка коэффициентов эллиптического оператора. Это означает априорное разбиение с гладкими границами, что практически не реализуемо па неструктурированных сетках для параллельных технологий типа "черный ящик".

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

Отправной точкой при построении декомпозиционного переобуславливателя является легко параллелезуемый аддитивный метод Шварца [199], представляющий собой просто блочпо-диагональную матрицу в случае минимального перекрывания подобластей. Диагональные блоки в ней — переобуславливатели к соответствующим диагональным блокам матрицы жесткости. Для оператора диффузии число обусловленности к таким образом переобусловленной матрицы жесткости зависит от диаметра подобластей 1/Н2), ширины перекрывания 1/5) и скачка коэффициента диффузии maxp/min/э). Если снабдить блочно-диагональный переобуславливатель коррекцией на каркасном подпространстве, оценка числа обусловленности улучшится то С(р) (l + у)2 [92], а в случае гладкого коэффициента р ее можно довести до С (l + у) [118, 238].

Здесь и далее C(z) обозначает константу, зависящую от z, но не зависящую от остальных параметров. Кроме этого, в дальнейшем с, С (с индексами и без индексов) будут обозначать некоторые положительные константы, не зависящие от параметров задачи, а выражение А ~ В будет обозначать спектральную эквивалентность матриц А и В с константами, не зависящими от размеров матриц, или приближенное равенство или пропорциональность величин А и В.

Отметим, что зависимость к от р не является критической, если коэффициент р(х) гладкий в больших подобластях и применяется метод сопряженных градиентов (СГ). Как показано в [150], при сильных скачках число итераций СГ лишь незначительно возрастает на величину, зависящую от геометрической структуры скачков р(х). Однако этот результат не применим для других итерационных методов, а также для гетерогенного коэффициента р(х).

В настоящей работе предлагается метод на основе аддитивного метода Шварца, впервые обеспечивающий теоретическую независимость к от гетерогенного коэффициента диффузии р(х) в зоне перекрывания и оценку к < С (l + j) (Теорема 1.3). В доказательстве теоремы используется гладкость р(х) в подобластях вне зоны перекрывания, хотя это не нужно на практике. Важным практическим достоинством метода является его независимость от гетерогенности р(х) в зоне перекрывания, что позволяет не координировать разрезание на подобласти с линией скачка. Суть метода — специфическое каркасное подпространство, получаемое агрегированием неизвестных внутри подобластей. Агрегирование неизвестных — весьма востребованная технология при разработке численных методов [86, 164, 253]. Предлагаемый метод, как и метод [164], использует негладкие базисные функции каркасного простраснтва, что отличает его от работ [86, 253]. В отличие от [164], наше каркасное пространство является гораздо большим, поскольку оно включает в себя базисные функции на мелкой сетке, которые не равны нулю в зоне перекрывания. С одной стороны, это делает каркасную (агрегированную) задачу более сложной для решения, с другой стороны, это позволяет исключить зависимость к, от р(х).

Параллелизация решения агрегированной задачи зависит от конкретного приложения. В случае двумерных сеток, уменьшение размерности агрегированной матрицы по отношению к исходной размерности настолько существенно (N —► iV1/2), что позволяет использовать прямой метод (факторизации) на каждом процессоре [78, 79]. В случае трехмерных сеток каркасная задача остается достаточно большой (понижение размерности N —► №/3), поэтому ее приближенное решение осуществляется с помощью блочного метода верхней последовательной релаксации BSOR или его симметризованного аналога. Эффективность подобного решения мотивируется тем, что арифметическая цена одной итерации BSOR существенно ниже умножения на исходную матрицу жесткости, а жесткость агрегированной матрицы существенно ниже жесткости исходной матрицы, что порождает высокую скорость сходимости BSOR. Легкая параллелезуемость метода является его привлекательной чертой, также как простота его реализации на неструктурированных сетках [187,188, 184]. По сути дела, метод является методом "серого ящика", поскольку помимо элементов исходной матрицы жесткости, требует только разбиение степеней свободы на непересекающиеся подмножества, обеспечиваемые, например, любым алгоритмом разбиения графа. Отметим, что по формальным признакам предлагаемый алгоритм относится к гибридным методам [194, 238], поскольку его аддитивная структура на уровне подобластей сочетается с мультипликативной процедурой коррекции в каркасном пространстве.

Доступные методы решения неструктурированных систем позволяют использовать полностью неструктурированные симплициальные сетки в адаптивных технологиях приближенного решения краевых задач, которые, в свою очередь, открывают путь к анизотропной адаптации, рассмотренной в главе второй. Как уже отмечалось выше, именно возможность анизотропной адаптации расширяет границы современных адаптивных технологий, поскольку позволяет эффективно аппроксимировать анизотропные решения. В ряде работ [111, 110, 112, 226], появившихся на рубеже 90-х годов, было показано, что треугольники с тупыми и острыми углами, вытянутые вдоль направления минимальной второй производной некоторой функции, могут быть наилучшими элементами для минимизации ошибки интерполяции. При этом степень вытянутости треугольников и их малости зависит от матрицы вторых производных (гессиана) функции. На практике управление вытянутостыо элементов задавалось метрикой на основе гессиана функции или его приближения [90, 117, 77, 213, 134]. Предпринимались попытки обобщений на адаптацию к вектор-функциям [116]. Однако первое теоретическое обоснование анизотропной адаптации появилось в работах [11, 53], где были доказаны асимптотические аппроксимационные свойства оптимальных симплициальных сеток и их приближений, квази-оптимальных сеток, которые можно строить на практике.

Технология адаптивного приближенного решения на основе восстановления гессиана сеточного решения породила ряд теоретических вопросов. На некоторые из них даны ответы во второй главе в виде теоретических результатов автора. Так, доказано существование оптимальных сеток при разумных предположениях на непрерывность ошибки интерполирования от координат сеточных узлов и не увеличении этой ошибки при иерархическом измельчении треугольников [12] (Теорема 2.3). Выведены асимптотические оценки ошибки в норме Ь^ для оптимальных сеток, в двумерном [11] и трехмерном [53] случаях, предложен их объединенный анализ [10] (Теорема 2.4), а в Теореме 2.7 приведены обобщения (совместно с А.Агузалом) на случай норм Ьр, 0 < р < +оо [10]. Рассмотрены свойства приближений оптимальных сеток, которые являются квазиравномерными в метрике, порожденной гессианом функции, и доказана Теорема 2.8 о том, что такие приближения являются квази-оптимальными сетками, ошибка иитерпо-ляции на которых асимптотически имеет такую же зависимость от числа симплексов сетки, что и на оптимальных сетках [11, 53].

Восстановление гессиана сеточной функции — важная часть адаптивной технологии. По сути дела, необходимо вычислить матрицу вторых производных у сеточного решения, что на практике осуществляется в рамках конечно-элементных аппроксимаций (за счет перебрасывания одной производной на пробную функцию в формуле Грина). Подобное восстановление не обеспечивает поточечную сходимость восстановленного гессиана к дифференциальному для функций с особенностями, поэтому предположение о поточечной аппроксимации, используемое в излагаемой теории, может не выполняться. С учетом этого автором был предложен новый метод восстановления [54] гессиана у функций с особенностями, который обеспечивает локальную аппроксимацию компонент дифференциального гессиана в Ь^ норме, что имеет смысл для функций из И^2. Здесь и далее обозначает соболевское пространство вещественных функций, у которых р-е степени производных порядка не выше q интегрируемы в рассматриваемой области. Отметим, что решения большинства эллиптических краевых задач принадлежат пространству И^2 в подобластях, где коэффициенты дифференциального оператора — гладкие функции [151]. Это открывает путь к локальной аппроксимации гессиана сеточных функций с особенностями. Доказанная автором Теорема 2.15 о локальной сходимости восстановленного гессиана к дифференциальному — первый теоретический результат для функций с особенностями и анизотропных сеток.

Генерация квази-равномерных сеток в заданной непрерывной метрике — ключевой элемент рассматриваемой адаптивной технологии [90, 117, 77, 213, 134]. В работе предлагается совместный алгоритм автора и К.Н.Липникова генерации тетраэдральных сеток, являющийся обобщением двумерного алгоритма [90] и представляющий собой последовательность локальных модификаций сетки для достижения ее квази-равпомерности.

Нужно отметить, что теоретического обоснования, почему сходится адаптивный алгоритм построения квази-оптимальных сеток для функций с особенностями, нет. Однако в Теореме 2.18 автору удалось построить оценки, показывающие возможность уменьшения нормы ошибки при интерполяции гладких функций [11, 53]. Вопрос теоретического обоснования в общем случае не решен, так же как и вопрос апостериорных оценок ошибки в рамках данной технологии, хотя для случая изотропного сгущения предложен конструктивный метод построения таких оценок [63], а для некоторых краевых задач и некоторых типов анизотропных сеток предложены оценки [174]. Альтернативная адаптивная технология [260, 59, 161], основанная на локальном измельчении/огрублении сетки на основе выравнивания апостериорно оцененной ошибки, естественным образом интегрирует оценку ошибки, однако использует иерахические локально измельчающиеся сетки, элементы которых наследуют форму начальной сетки.

Общий адаптивный алгоритм весьма технологичен, поскольку его структура предполагает выполнение четырех совершенно независящих друг от друга шагов:

1. шаг инициализации: построение начальной сетки;

2. нахождение сеточного решения (непрерывной кусочно-линейной функции);

3. восстановление сеточного гессиана;

4. генерация сеток, квази-равномерных в заданной метрике.

Ключевым свойством алгоритма является полная автономность от пользователя шагов 3 и 4, реализуемых по принципу "черных ящиков", а также минимизация обмена данными между сетками (используются только данные о сетке и сеточном решении). Удобство использования метрики заключается в том, что ее легко модифицировать, тем самым управляя процессом адаптации и свойствами симплициальной сетки. Например, можно фокусировать измельчение сетки в зонах повышенного интереса и ослаблять концентрацию узлов сетки в зонах пониженного интереса. Аналогичные цели типичны для целевых адаптивных технологий [219]. Кроме того, можно управлять степенью вытянутости симплексов, если используемый метод аппроксимации не может использовать анизотропные сетки. В Теореме 2.19 автором показано, как влияет модификация метрики на ошибку интерполяции, через сравнение ошибок на квази-оптимальных сетках и сеток с управляемой адаптацией [13, 190] ; там же рассмотрено несколько примеров управления.

Восстановление гессиана сеточной функции нашло свое применение и в построении поверхности более высокого порядка на основе кусочно-линейных поверхностей. Актуальность этой проблемы вызвана тем, что современные системы САПР порождают границы расчетной области в виде треугольных сеток в пространстве. Недостаточное разрешение поверхности границы зачастую ухудшает работу адаптивных методов, поскольку дискретная граница сама по себе порождает дополнительную ошибку аппроксимации, не уменьшаемую за счет адаптации внутри области. Точность аппроксимации кусочно-гладкой поверхности можно существенно повысить за счет восстановления поверхности более высокого порядка, чем заданная кусочно-линейная дискретизация. Существует несколько методов такого восстановления [142, 200, 203, 204]: в одних [200, 204] параметризация поверхности (а также другие характеристики поверхности) вычисляется из производных функций, задающих параметризацию; в других [142, 200] дискретная кусочно-линейная поверхность аппроксимируется кусочно-квадратичной поверхностью методом наименьших квадратов. Предлагаемый в диссертации метод использует технологию вычисления приближенного гессиана кусочно-квадратичной функции, задающей восстанавливаемую поверхность. Гессиан вычисляется на основе обобщенной формулировки, по аналогии с методом конечных элементов. Предлагаемый метод точен для квадратичных поверхностей и позволяет повысить порядок аппроксимации кусочно-гладкой границы, если ее куски представимы функциями с непрерывными третьими производными. В Теореме 2.22 доказаны оценки аппроксимации для кусочно-гладких границ общего вида.

Приближенное решение трехмерных краевых задач даже на адаптивных сетках требует весьма большого числа симплексов, поэтому вопросы параллельной адаптации наиболее актуальны для трехмерных задач. Генерация трехмерной сетки, квази-равномерной в заданной метрике, является весьма трудоемкой процедурой, и может занимать больше времени, чем получение сеточного решения. Поэтому ускорение приближенного адаптивного решения требует параллелизации как решения сеточных систем, рассмотренной в главе первой, так и генерации трехмерных сеток. Предлагаемый метод параллельной генерации основан па разрезании сетки на подобласти-слои с чередующимися направлениями нарезки, независимой генерации сеток в слоях с замороженными границами, и склейке подсеток в глобальную сетку [187, 189, 191]. Численные расчеты для некоторых модельных задач математической физики и гидродинамики [187, 188] показали параллельную эффективность предлагаемой технологии адаптивного решения, выявив при этом интересную особенность параллельной генерации: если сетки почти установились в ходе адаптации, а число локальных модификаций сетки достаточно велико, то наблюдается сверхлинейное параллельное ускорение процесса генерации сетки [187, 189, 191].

Разбиение расчетной области на подобласти и использование в них конформных сеток, не согласованных друг с другом на границах подобластей, порождает пестыкую-щиеся сетки, впервые проанализированные в [70]. Технологические преимущества таких сеток были перечислены выше. Анализ аппроксимациопных свойств пестыкующихся сеток в рамках метода мортарных конечных элементов и/или макро-гибридных постановок представлен в работах [52, 67, 68, 70, 169, 73] и сводится к утверждению, что при корректной трактовке условий сшивки на интерфейсе асимптотические свойства энергетической нормы ошибки аналогичны результатам для конформных сеток.

В этой связи уместно подчеркнуть исключительную важность корректности условий сшивки решения на интерфейсе между подобластями, которые неразрывно связаны со свойствами операторов Пуанкаре-Стеклова. На функциональном уровне условия сшивки и свойства операторов Пуанкаре-Стеклова в приложении к методам декомпозиции области были впервые рассмотрены в монографиях [1, 31], а также в работах [2, 50, 51]. Общность полученных результатов породила ряд композиционных методов для задач с различными операторами в неперекрывающихся подобластях [31, 220]. Однако эта методология не переносится на случай нестыкующихся сеток и на случай многих подобластей с произвольными скачками коэффициентов. Эффективные для этих случаев алгоритмы базируются на специальных интерфейсных переобуславливателях, строящихся на сеточном, а не дифференциальном уровне.

Системы сеточных уравнений, порождаемые аппроксимациями на нестыкующихся сетках, являются седловыми. Начиная с середины девяностых годов, эффективному решению этих систем было посвящено много работ [39, 69, 44, 47, 46, 178, 179, 239, 160, 128, 127, 122, 181, 256, 80, 148, 271, 268]. Эти методы можно условно разделить на три группы: многосеточные методы, использующие иерархию сеток как в подобластях, так и на интерфейсах между подобластями [80, 148, 271, 268], декомпозиционные методы, строящие интерфейсные переобуславливатели на основе приближенного обращения дополнения Шура на разрежающихся сетках [44, 46, 178, 179, 163, 160, 128, 127, 122], а также декомпозиционные интерфейсные переобуславливатели [47, 239, 181, 256]. Достоинством методов первой группы является независимость скорости сходимости от шага сетки и асимптотическая оптимальность арифметической работы, а недостатком — требование иерархичности сеток; достоинством второй группы также является независимость скорости сходимости от шага сетки, причем технология может как опираться [44, 46, 178, 179, 160, 128, 127, 49], так и не опираться [122] на иерархичность сетки, а недостатком - зависимость арифметической работы от структуры сгущения исходной сетки и достаточно трудоемкое применение интерфейсного переобуславливателя; достоинством третьей группы [239, 181, 256] является удобство реализации (технологичность) и параллелизации, а недостатком — возможная незначительная зависимость скорости сходимости от шага сетки. В третьей главе автором рассмотрены методы второй и третьей группы с точки зрения их применения в адаптивном решении краевых задач. В Лемме 3.5 и Теореме 3.7 обоснована независимость скорости сходимости итерационных алгоритмов от шага сетки (методы второй группы), а также числа подобластей и величины коэффициентов диффузии и реакции для уравнений реакции-диффузии (методы обеих групп). Построены параллельные адаптивные алгоритмы и на численных примерах проиллюстрированы их основные свойства.

Аппроксимации на нестыкующихся сетках являются неконформными, поскольку используемые конечно-элементные пространства порождают разрывные на интерфейсах решения, которые не принадлежат стандартному пространству обобщенных решений И^П). Специфические варианты подобных аппроксимаций удобно применять для декомпозиционных алгоритмов и на стыкующихся на интерфейсе сетках. Действительно, в качестве условий сшивки на интерфейсе между стыкующимися сетками можно потребовать поузловую непрерывность сеточного решения и его потока на интерфейсе [131] и сформулировать и обосновать эффективные алгоритмы итерационного решения получаемых седловых систем [252]. С учетом конформности сетки можно показать, что полученное решение удовлетворяет стандартной конечно-элементной системе на стыкующихся сетках. Для некоторых случаев разбиения области на подобласти удается использовать нормировку сеточных функций в пространстве следов [205, 207] на интерфейсах в качестве спектрально эквивалентного переобуславливателя, при этом эффективность его применения достигается технологией мозаично-скелетонной аппроксимации [249, 251]. Основная идея состоит в аппроксимации плотной матрицы, возникающей из нормализации следа, некоторой разреженной матрицей, которую легко умножать на вектор. Для двумерного метода подструктур такая попытка предпринималась в [171]. В работе [252] рассматривается интерфейсное переобуславливание дуального дополнения по Шуру для множителей Лагранжа, что оказывается гораздо проще даже в трехмерном случае. Предлагаемый переобуславливатель может быть реализован как "серый ящик", требующий на входе коэффициенты задачи и интерфейсную сетку, поэтому напоминает быстрые муль-типольные методы [66].

Другим источником неконформности метода даже на конформных сетках являются конечно-элементные пространства, базисные функции которых не принадлежат пространству И^П). Практическая значимость простейшего из них, Крузэ-Равиара [108], определяется тем обстоятельством, что порождаемые конечно-элементные матрицы имеют тот же вид, что и статически конденсированные матрицы смешанных конечных элементов Равиара-Тома наименьшего порядка [225], обеспечивающие локально консервативные аппроксимации эллиптических уравнений на неструктурированных сетках и одновременное вычисление скалярного решения и его потока. Параллельная реализация [102] многосеточного метода для конечно-элементных аппроксимаций Крузэ-Равиара на локально сгущающихся иерархических симплициальных сетках позволила построить эффективную параллельную технологию [103] решения смешанных конечно-элементных систем:

1. на входе дается линейная смешанная конечно-элементная система в дисассемблированном (по ячейкам сетки) виде, а также координаты вершин симплексов и граф связности вершин;

4. итерационно решается разреженная конденсированная система для множителей Лагранжа;

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

Ключевым наблюдением является то, что все этапы, кроме четвертого, арифметически дешевы и легко параллелезуемы, в силу локальности действий. Параллельное решение линейной системы, как наиболее трудоемкий этап всей цепочки, осуществляется быстрос-ходящимся многосеточным методом [102, 103, 16], реализация которого позволила сформировать всю технологическую цепочку.

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

Алгоритмы четвертой главы объединяет применение и тестирование в приложениях на прямоугольных сетках, хотя использование ни одного из методов не ограничено классом прямоугольных сеток. Рассматриваются две основные задачи гидродинамики: уравнения Навье-Стокса и задача многофазной фильтрации. Являясь нелинейными нестационарными уравнениями, они допускают ряд технологий получения решения [216, 132, 156, 243, 198, 145, 3, 26, 4]. Здесь мы ограничимся двумя наиболее востребованными методами. Во-первых, это полностью неявные аппроксимации по времени, обеспечивающие абсолютную устойчивость к шагу по времени и порождающие на каждом временном шаге систему нелинейных уравнений, которую необходимо решить. Во-вторых, это проекционный алгоритм для уравнений Навье-Стокса, сводящийся к решению одной или нескольких линейных систем, с нежесткими ограничениями на шаг по времени [198, 154, 248]. Преимуществом первого метода является возможность управлять шагом по времени на основе внешних данных (правой части, граничных условий), а не внутренней динамики системы, что весьма привлекательно в инженерных приложениях. Преимущества второго метода — в упрощении решаемых систем уравнений и существенном ускорении их решения на каждом шаге, что обеспечивает этим методам наибольшую скорость расчетов [248].

Проекционные алгоритмы были предложены в работах [101, 242, 42] и их смысл сводился к двухшаговой процедуре: с помощью уравнения моментов предсказать поле скоростей на очередном шаге по времени и спроектировать это поле в пространство соле-иоидальных функций. Предсказание возможно за счет явной аппроксимации по времени уравнения моментов [101] (при этом налагаются существенные ограничения на временной шаг), за счет метода характеристик [216] (параллельная реализация которого очень сложна [45]), за счет полунеявной аппроксимации по времени уравнения моментов [198, 170], чье единственное отличие от неявной аппроксимации — линеаризация нелинейного конвективного члена вида (и'с+1 • У)и'с+1 —> (и*2 • У)и'с+1. Последний подход предполагает решение линеаризованного уравнения моментов, представляющего собой в трехмерном случае три сингулярно возмущенных уравнения реакции-конвекции-диффузии для каждой из компонент скорости. Именно это приложение потребовало разработки эффективного параллельного метода решения сингулярно возмущенного уравнения конвекции-диффузии.

Сингулярно возмущенные операторы предоставляют дополнительные возможности для эффективной параллелизации решателей на основе методов декомпозиции с перекрывающимися подобластями. В применении к параболическим уравнениям наличие сингулярно возмущенного члена нулевого порядка привело к ряду родственных методов [176, 177, 74], основанных на экспоненциальном убывании сеточной функции Грина.

Предлагаемый алгоритм базируется на важном свойстве сингулярно возмущенного уравнения конвекции-диффузии: возмущения решения, вызванные локальным возмущением краевого условия, распространяются анизотропно. В случае постоянного вектора переноса точечное возмущение распространяется только вниз по потоку, быстро убывая в поперечном направлении и вверх по потоку. Это свойство может быть переформулировано в терминах анизотропного поведения сеточной функции Грина [165, 211, 137], которое может быть использовано при построении как декомпозиционных [222, 137], так и многосеточных методов [35, 212].

На основе этих свойств оператора конвекции-диффузии конструируется параллельный двухуровневый метод Шварца [137, 141], быстро сходящийся даже для малых налеганий между подобластями. На первом уровне область разбивается на налегающие полосы (слои), перпендикулярные доминирующему направлению переноса, и определяются итерации Шварца в направлении вниз по потоку. На втором уровне каждая подобласть-полоса разбивается на налегающие подобласти-квадраты и задаются несколько итераций Шварца между квадратами с четными и нечетными индексами. Естественный параллелизм методу дает именно второй уровень, поскольку каждая подобласть-квадрат может обрабатываться своим процессором. Если вектор переноса сильно отклоняется от заданного направления (например, вследствие завихрений), разбиение первого уровня должно быть адаптировано к поведению вектора.

Рассмотренный метод является первым итерационным методом для оператора конвекции-диффузии с плоско-параллельпым полем переноса, чья скорость сходимости теоретически ограничена сверху константой, не зависящей от малости коэффициента диффузии е и шага сетки к, если £ < к, см. Теорему 4.1, впервые доказанную в [136]. Метод обобщен автором на случай полей переноса с локальными завихрепностями, при этом в Теореме 4.6 автором доказано, что итерационная скорость сходимости остается не зависящей от малости е и к при разумных предположениях на разбиение на подобласти. Арифметическая цена одной итерации этого метода (как и его инициализации) пропорциональна числу неизвестных в системе, поскольку представляет собой набор линейных систем с матрицами ограниченного порядка (в двумерном случае) или ленточными матрицами с малой шириной ленты (в трехмерном случае). Более того, метод легко параллелезуем. Отметим, что предложенный метод, будучи первым проанализированным декомпозиционным алгоритмом, универсальным по отношению к сеточному числу Пекле, получил обобщения на нелинейные задачи [75, 76]. Добавление в уравнение сингулярно возмущенного члена нулевого порядка (реакции) только улучшает свойства метода, поскольку позволяет избегать адаптивного разбиения для поля переноса с завихрениями, сохраняя высокую скорость сходимости [138]. Среди методов декомпозиции на неперекрывающиеся области отметим метод [144, 48] со скоростью сходимости, не зависящей от малости коэффициента диффузии е.

Параллельные приложения предлагаемого метода декомпозиции в рамках проекционного алгоритма весьма разнообразны: это решение нестационарных уравнений Навье-Стокса [198, 248], решение уравнения тепло-массопереноса в простейшем химическом реакторе [166], моделирования взаимодействия течения вязкой несжимаемой жидкости и упругой стенки канала [215, 233]. Общая черта всех проекционных алгоритмов — решение на каждом временном шаге сеточного аналога уравнения Пуассона для поправки к давлению, которая используется при проектировании предсказанного поля скоростей на множество функций с заданной дивергенцией. Решения последовательности этих жестких задач должны быть найдены с высокой точностью, причем физически мотивированное начальное приближение является нулевым вектором. Выбор начального вектора, более близкого к решению очередной системы, позволяет существенно сократить количество итераций и расчетное время.

Технология выбора начального вектора для систем с фиксированной симметричной положительно определенной матрицей и последовательностью правых частей рассматривалась в [95, 129], а для систем с фиксированной несимметричной матрицей — в [247]. Основной идеей технологии является использование объединенного крыловского подпространства для всех систем. Особенности аппроксимации уравнений Навье-Стокса [140, 247] могут породить системы с несимметричными матрицами, поэтому мы ограничимся рассмотрением несимметричных систем. Начальное приближение в этих методах, — решение спроектированной на крыловское подпространство очередной системы, — оказывается очень хорошим при условии достаточной информативности подпространства. Накапливание крыловских векторов — необходимое условие реализации технологии — проводит к повышенным запросам на память компьютера, что удовлетворяется параллельными компьютерами с распределенной памятью. В некоторых приложениях [166] проекционный алгоритм порождает последовательность систем с меняющимися матрицами. В таких случаях накапливание объединенного крыловского подпространства заменяется аккумуляцией независимых пространств (для каждой системы) и последовательным проектированием текущего начального приближения на каждое из этих пространств [245]. Эффективность такого проектирования зависит от того, насколько сильно отличаются друг от друга спектральные базисы матриц из последовательности систем. Отметим в этой связи методику выбора векторов из последовательности крыловских подпространств для ускорения итерационной скорости сходимости [149, 227].

Решение последовательности нелинейных систем уравнений, возникающей из полностью неявной аппроксимации нестационарных уравнений, также может быть ускорено за счет выбора более близкого начального приближения. Более того, квадратичная скорость сходимости метода Ньютона повышает важность такого выбора. Физически мотивированный выбор начального приближения — использование решения с предыдущего временного шага. Методы построения начального приближения на основе нескольких предыдущих решений рассматривались в [133].

В настоящей работе рассматривается новый алгоритм INB-POD [246, 245], базирующийся на собственном ортогональном разложении (proper orthogonal décomposition POD), которое порождает оптимальные аппроксимации малых размерностей для наборов данных. POD генерирует ортонормированный базис для представления данных с помощью метода наименьших квадратов [223, 224]. Используя галеркинские проекцию текущего нелинейного оператора на базис POD, можно построить текущую нелинейную модель гораздо меньшей размерности. Эти модели дают лучшее начальное приближение для метода Ныотона на следующем временном шаге. Платой за лучшее приближение является дополнительная нелинейная система, которую необходимо решать на каждом временном шаге. Однако арифметическая цена этого решения оказывается меньшей, чем решение исходной системы (в силу малой размерности), а уменьшение ньютоновских итераций настолько значительно, что общее время решения уменьшается в несколько раз.

Важным преимуществом алгоритма INB-POD в параллельных вычислениях (стандартных или распределенных) является то, что он не строго детерминирован: неявная схема может считаться как без, так и с POD ускорением, а также как с новой, так и устаревшей упрощенной моделью, причем POD базис может вычисляться на другом процессоре. Технология асинхронных обменов позволяет использовать обновленную упрощенную модель только когда соответствующий базис будет доступен, т.е. вычислен где-то и получен процессором, считающим неявную схему. Со своей стороны, этот процессор по окончании каждого временного шага может асинхронно отсылать решение POD генератору. Такая технология идеально подходит для распределенных вычислений (meta-computing) [246], осуществляемых на удаленных вычислительных системах с медленной и неустойчивой связью. Предлагаемая технология тестировалась на полностью неявном решении двумерного уравнения Навье-Стокса в формулировке для функции линий тока (задача о каверне).

В отличие от уравнений Навье-Стокса, уравнения многофазной фильтрации [3, 26, 265] характеризуются сложными нелинейными зависимостями, гетерогенностью и возможным вырождением коэффициентов уравнения, а также наличием дополнительных уравнений, описывающих течения в скважинах. Не касаясь теории дифференциальных уравнений, перейдем к обсуждению методов их приближенного решения. С повышением производительности расчетов, полностью неявные аппроксимации по времени сгали общепринятыми в индустриальном сообществе. Как правило, якобиан порожденной нелинейной системы известен аналитически, поэтому неточный метод Ньютона интегрирован в программу моделирования [184, 267, 124, 26]. С практической точки зрения, повышение итерационной скорости сходимости решения линейных систем с якобианом является очень важным в деле повышения эффективности всего моделирования. В случае многофазных многокомпонентных течений, линейная система содержит уравнения, порожденные различными процессами, происходящими в разных масштабах, поэтому различные блоки якобиевой матрицы могут иметь различные математические свойства [173]. В силу этого, прямой перенос современных технологий (многосеточные методы, методы декомпозиции), разработанных для эффективного решения эллиптических уравнений, невозможен. Разумеется, для упрощенных моделей, требующих решения сеточных эллиптических задач, миогие технологии были успешно перенесены [262, 263, 167, 266], однако решение многокомпонентных систем требует дополнительных усилий.

Одним из наиболее эффективных подходов оказался метод алгебраического выделения уравнения для давления [173, 261, 124, 185, 184], который преобразует исходную систему к специальному блочному виду, стремясь максимально отщепить зависимость переменной "давление" от переменных "концентрации". Такой подход имеет два потенциальных преимущества: во-первых, можно строить хороший переобуславливатель только для самых жестких блоков, переобуславливая другие блоки вычислительно дешевыми технологиями. Во-вторых, этот подход не зависит от структуры сетки, поскольку пере-обуславливатели для блоков могут опираться только на их алгебраическую структуру.

Второй обсуждаемый подход — многосеточный метод с ускоренным разгрублением БОУЮ [267, 184], который использует иерархичность трехмерной прямоугольной сетки. Несмотря па структурное подобие БСА/Ю и стандартного У-цикла геометрического многосеточного метода Мв [40, 158], они опираются на разные итерационные механизмы. Трехмерный метод МО использует последовательные коррекции на более грубых уровнях, в то время как БСМС — переобуславливатель без пост-сглаживания с коррекцией (аналогичной методу агрегирования из главы 1) в подпространстве двумерных функций, где используется усиленный (со многими сглаживаниями, в силу малости задач) двумерный метод Мв для решения агрегированной подзадачи. Именно это обстоятельство позволяет использовать в БСМС простейшие геометрические операторы интерполяции и сужения даже для приложений с сильно меняющимися коэффициентами, что является неприемлемым для классического метода трехмерного метода МО [115]. Несмотря на разные механизмы, лежащие в основе используемых методов (возможность отщепления давления и эффективность вертикального агрегирования), численные эксперименты показали эффективность обоих методов решения возникающих линейных систем.

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

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

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

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

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

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

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

В разделе 2.1 представлен теоретический анализ этих сеток: введены основные понятия, доказано существование оптимальных сеток в двумерном случае, для норм ошибки показаны оценки снизу на оптимальных сетках как в норме Loo > так и в норме Lp, р € (0, +оо], и оценки сверху на квази-оптимальных сетках в норме L^.

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

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

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

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

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

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

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

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

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

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

В разделе 4.1 формулируется и анализируется параллельный двухуровневый метод Шварца в применении к сингулярно возмущенному уравнению конвекции-диффузии.

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

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

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

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

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

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

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

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

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

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

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

Основные результаты диссертации докладывались автором: на всероссийской школе-конференции (Казань, 1999), всероссийских конференциях по построению расчетных сеток (Москва, 2002,2004), международных конференциях по методам декомпозиции области DDM (Чиба, 1999), параллельной вычислительной гидродинамике ParCFD (Москва, 2003), геофизическим наукам SIAM (Остин, 2003, Авиньон, 2005), Европейских конференциях по вычислительной математике ENUMATH (Париж, 1995, Юваскюла, 1999, Ис-кия, 2001), Европейской конференции по вычислительным методам ECCOMAS (Париж, 1996), Франко-Русско-Итало-Узбекских симпозиумах по численному анализу и приложениям (Ташкент, 1995, Марсель, 1997), международной конференции по анализу, вычислениям и применениям дифференциальных и интегральных уравнений (Штуттгарт, 1996), международной конференции GAMM по параллельным многосеточным методам (С.-Вольфганг, 1996), международной конференции "Domain decomposition and multifield theory" (Обервольфах, 1998), международных симпозиумах "Finite element rodeo" (Хьюстон, 1999, Колледж Стэйшн, 2002) международной конференции "50 лет сопряженным градиентам"(Юваскюла, 2002), на научно - исследовательских семинарах Института вычислительной математики РАН, Вычислительного центра РАН, Института математического моделирования РАН, Университетов Париж 6, Париж 13, Лион 1, Ренн, Гейдель-берг, Мюнхен, Аугсбург, Юваскюла, Наймеген, Остин, Хьюстон, INRIA, Institut Francais du Pétrol, ExxonMobil Upstream Research С.

По теме диссертации опубликовано 42 работы, из них 21 в рецензируемых журналах, 18 в материалах конференций, 3 в научных изданиях. Вклад автора в совместные работы заключался: в формировании постановки проблемы [121, 11, 53, 185, 54, 187, 184, 12, 102, 189, 10, 123, 13, 15, 246, 103, 245, 122, 252, 190, 191, 188, 14, 16], предложении идеи решения [121, 11, 53, 185, 137, 140, 54, 187, 184, 12, 102, 189, 10, 123, 13, 15, 246, 103, 245, 122, 252, 49, 190, 191, 161, 188, 78, 8, 14, 16, 79, 139, 182], теоретическом обосновании [121,11,185,187,184,12,15, 246,103,122,181,14], совместном теоретическом обосновании [160, 53, 137, 140, 54, 102, 189, 10, 123, 13, 245, 138, 252, 190, 191, 78, 8, 79, 139, 182], технической реализации [137, 140, 15, 245, 141, 138, 181, 139], совместной технической реализации [121, 160, 11, 53, 185, 187, 184, 102, 123, 246, 103, 163, 183, 122, 128, 127, 252, 49, 189, 13, 190, 191, 161, 188, 78, 14, 16, 79, 162, 182], постановке экспериментов [121, 160, 11, 53, 137, 185, 140, 187, 184, 102, 189, 123, 13, 15, 246, 103, 245, 163, 183, 122, 128, 141, 127, 138, 181, 252, 49, 190,191, 188, 78, 14, 79, 139, 162, 182]. Содержание раздела 1.3 базируется на экспериментах автора с неопубликованным методом, предложенным и реализованным С.А.Горейновым, исходя из постановки проблемы автором.

Автор искренне признателен всем соавторам и коллегам за сотрудничество. Автор сердечно благодарит Ю.А.Кузнецова за многолетнюю поддержку, и С.В.Непомнящих за плодотворное обсуждение многих научных вопросов.

Заключение диссертации по теме "Вычислительная математика"

Выводы главы 4

Перечислим основные результаты четвертой главы.

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

2. рассмотрены численные параллельные трехмерные приложения разработанной технологии (нестационарные уравнения Навье-Стокса, уравнения тепло-массопереноса в простейшем химическом реакторе);

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

Заключение

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

Список источников диссертации и автореферата по математике, доктора физико-математических наук, Василевский, Юрий Викторович, Москва

1. Агошков В.И., Лебедев В.И. Операторы Пуанкаре-Стеклова и их приложение в анализе. -М.: ОВМ АН СССР, 1983.

2. Агошков В.И., Лебедев В.И. Операторы Пуанкаре-Стеклова и методы разделения области в вариационных задачах. // Выч.процессы и системы. -М.: Наука, 1985, с.173-227.

3. Азиз X., Сеттари Э. Математическое моделирование пластовых систем. -М.: Недра, 1982.

4. Басниев К.С., Дмитриев Н.М., Каневская Р.Д., Максимов В.М. Подземная гидромеханика. М.-Ижевск: Институт компьютерных исследований, 2005.

5. Бахвалов Н.С. О сходимости одного релаксационного метода для эллиптического оператора с естественными ограничениями.// ЖВМ и МФ, 1966, Т.6, с.101-135.

6. Бахвалов Н.С. Эффективные методы решения жестких многомерных многопараметрических задач. // ЖВМ и МФ, 1999, Т.39, с.2019-2049.

7. Булеев Н.И. Численный метод решения двумерных уравнений диффузии. // Атомная энергия. 1959, Т.6, N3, с.338-340.

8. Василевский Ю.В., Ильин В.П., Тыртышников Е.Е. Вычислительные технологии, // Современные проблемы вычислительной математики и математического моделирования, Т.1, М.: Наука, 2005, с.100-148.

9. Василевский Ю.В. Методы решения краевых задач с использованием нестыкую-щихся сеток. // Тр. Матем. центра им. Н.И.Лобачевского, Т.2. Итерационные методы решения линейных и нелинейных сеточных задач. Казань: УНИПРЕСС, 1999, с.94-121.

10. Василевский Ю.В., Агузал А. Объединенный ассимптотический анализ интерполяционных ошибок на оптимальных сетках. // ДАН, 2005, Т.405, N0.3, с. 1-4.

11. Василевский Ю.В., Липников К.Н. Адаптивный алгоритм построения квазиоптимальных сеток.// ЖВМ и МФ, 1999, Т.39, No.9, с.1532-1551.

12. Василевский Ю.В., Липников К.Н. Оптимальные триангуляции: существование, аппроксимация и двойное дифференцирование Р\ конечно-элементных функций. // ЖВМ и МФ, 2003, Т.43, No.6, с.866-874.

13. Василевский Ю.В., Липников К.Н. Оценки ошибки для управляемых адаптивных алгоритмов на основе восстановления гессиана. // ЖВМ и МФ, 2005, Т.45, No. 8, с.1424-1434.

14. Василевский Ю.В., Прокопенко A.B. Факторизация сеточного лапласиана на последовательностях сеток: экспериментальные оценки вычислительной работы.// Методы и технологии решения больших задач. -М.: ИВМ РАН, 2004, с.103-118.

15. Василевский Ю.В., Чугунов В.Н. О параллельном многосеточном решении консервативных неструктурированных конечно-элементных систем. // Методы и технологии решения больших задач. -М.: ИВМ РАН, 2004, с.91-102.

16. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. -М: Наука, 1984.

17. Воеводин В.В., Воеводин В.В. Параллельные вычисления. -Петербург: БХВ-Петербург, 2004.

18. Горейнов С.А. Мозаично-скелетонные аппроксимации матриц, порожденных асимптотически гладкими и осциляционными ядрами. // Матричные методы и вычисления, -М: ИВМ РАН, 1999, с.42-72.

19. Деммель Дж. Вычислительная линейная алгебра. -М.: Мир, 2001.

20. Джордж А., Лю Дж. Численное решение больших разреженных систем уравнений. -М.: Мир, 1984.ч 22. Дьяконов Е.Г. Минимизация вычислительной работы. Асимптотически оптимальные алгоритмы для эллиптических задач. -М.: Наука, 1989.

21. Иваненко С.А. Адаптивно-гармонические сетки. -М.: ВЦ РАН, 1997.

22. Ильин В.П. Методы неполной факторизации для решения алгебраических систем, -М.: Наука, 1995.

23. Ильяш Ю.И. Алгебраические многосеточные методы и их приложения к модельным задачам вычислительной гидродинамики. Дисс. к.ф.-м.н. -М.: ИВМ РАН, 1996.

24. Каневская Р.Д. Математическое моделирование гидродинамических процессов разработки месторождений углеводородов. -М.-Ижевск, Институт компьютерных исследований, 2003.

25. Кобельков Г.М. О численных методах решения уравнений Навье-Стокса в переменных скорость-давление. // Вычислительные процессы и системы, Вып.8, -М.: Наука, 1991, с.204-236.

26. Кузнецов Ю.А. Вычислительные методы в подпространствах. // Вычислительные процессы и системы, Вып.2, -М.: Наука, 1985, с.265-350.

27. Марчук Г.И., Лебедев В.И. Численные методы в теории переноса нейтронов. -М.: Атомиздат, 1981.

28. Лебедев В.И. Методы композиции. -М.: ОВМ АН СССР, 1986.

29. Марчук Г.И. Методы вычислительной математики. -М.: Наука, 1989.

30. Оганесян J1.A., Руховец Л.А. Вариационно-разностные методы решения эллиптических уравнений. -Ереван: Изд. АН Арм.ССР, 1979.

31. Ольшанский М.А. Лекции и упражнения по многосеточным методам. -М: Физмат-лит, 2005.

32. Ольшанский М.А. Анализ многосеточного метода для уравнений конвекции-диффузии с краевыми условиями Дирихле.// ЖВМ и МФ, 2004, Т.44, N8, с. 14501479.

33. Писсанецки С. Технология разреженных матриц. -М.: Мир, 1988.

34. Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. -М.: Наука, 1978.

35. Федоренко Р.П. Релаксационный метод решения разностных эллиптических уравнений.// ЖВМ и МФ, 1961, T.I, No.5, с.922-927.

36. Чижонков Е.В. Релаксационные методы решения седловых задач. -М.: ИВМ РАН, 2002.

37. Шайдуров В.В. Многосеточные методы конечных элементов, -М.: Наука, 1989.

38. Шишкин Г.И. Сеточные аппроксимации сингулярно возмущенных и параболических уравнений. -Екатеринбург: УО РАН, 1992.

39. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. -Новосибирск: Наука, 1967.

40. Abakumov A., Yeremin A., Kuznetsov Yu. Efficient fast direct method of solving Poisson's equation on a parallelepiped and its implementation in an array processor. // Sov.J. Numer.Anal. Math.Modelling, 1988, V.3, p.1-20.

41. Abdoulaev G., Achdou Y., Kuznetsov Yu., Prud'homme C. On a parallel implementation of the mortar element method.// Math. Model. Numer. Anal. (M2AN), 1999, V.33, p.245-259.

42. Abdoulaev G., Achdou Y., Hontand J., Kuznetsov Y., Pironneau 0., Prud'homme C. Nonmatching grids for fluids. // Proc. 10th Int. Conf. on DDM for PDE's. Cont. Math., 1998, V.218, p.3-22.

43. Achdou Y., Kuznetsov Yu., Pironneau O. Substructuring preconditioners for the Q^ mortar element method.// Numer.Math., 1995, V.71, p.419-449.

44. Achdou Y., Maday Y., Widlund O. Iterative substructuring preconditioners for mortar element methods in two dimensions.// SIAM J. Numer. Anal., 1999, V.36, No.2, p.551-580.

45. Achdou Y., Le Tallec P., Nataf F., Vidrascu M. A domain decomposition preconditioner for an advection-diffusion problem.// Сотр.Methods Appl.Mech.Engnr., 2000, V.184, p.145-170.

46. Achdou Y., Jaffre J., Svyatski D., Vassilevski Yu. Numerical simulation of unsteady 3D flows on anisotropic grids. // Transactions of French-Russian Liapounov Institute for Applied Mathematics and Computer Science, V.4, Moscow, 2003, p.107-124.

47. Agoshkov V., Lebedev V. Generalized Schwartz algorithm with variable parameters. // Sov. J. Numer. Anal. Math. Modelling, 1990, V.5, p.1-26.

48. Agoshkov V., Lebedev V. Variational algorithms of the domain decomposition method. // Sov. J. Numer. Anal. Math. Modelling, 1990, V.5, p.27-46.

49. Agouzal A., Thomas J.-P. Une méthode d'élement finis hybridesen décomposition de domaines. // RAIRO M2 AN, 1995, V.29, p.749-764.

50. Agouzal A., Lipnikov K., Vassilevski Yu. Adaptive generation of quasi-optimal tetrahedral meshes. // East-West J. Numer. Math., 1999, V.7, p.223-244.

51. Agouzal A., Vassilevski Yu. On a discrete Hessian recovery for Pj finite elements.// J. Numer. Math., 2002, V.10, p.1-12.

52. Amestoy P., Duff I., L'Excellent J., Li X. Analysis and comparison of two general sparse solvers for distributed memory computers. // ACM Transactions on Mathematical Software, 2001, V.27, p.388-421.

53. Apel T., Randrianarivony H. Stability of discretizations of the Stokes problem on anisotropic meshes. // Math. Comput. Simulât., 2003, V.61, p.437-447.

54. Apel T. Anisotropic finite elements: Local estimates and applications. -Stutgart: Teubner, 1999.

55. Axelsson O. Iterative solution methods. -N.Y.: Campridge Univ.Press, 1996.

56. Babuska I., Stroubolis T. The Finite Element method and its Reliability. -Oxford: Clarendon Press, 2001.

57. Bânsch E. Local mesh refinement in 2 and 3 dimensions. // IMPACT of Computing in Science and Engrg., 1991, V.3, p.181-191.

58. Banegas A. Fast Poisson solvers for problems with sparsity.// Math.Comp., 1978, V.32, p.441-446.

59. Bank R., Rose D. Marching algorithms for elliptic boundary value problems, I: The constant coefficient case.// SIAM J. Numer. Anal., 1977, V.14, p.792-829.

60. Bank R., Xu J. Asymptotically Exact A Posteriori Error Estimators, Part II: General Unstructured Grids. // SIAM J. on Numer. Anal., 2004, V.41, No.6, p.2313-2332.

61. Barrett R., Berry M., Chan T., Demmel J., Donato J., Dongarra J., Eijkhout V., Pozo R., Romine C., Van der Vorst H. Templates for the solution of linear systems: building blocks for iterative methods. -Philadelphia, PA: SIAM, 1994.

62. Barth T. Aspects of unstructured grids and finite volume solvers for the Euler and Navier-Stokes equations. AGARD Report 787, 1992.

63. Beatson R., Greengard L. A short course on fast multipole methods. // Wavelets, Multilevel Methods and Elliptic PDEs (Eds. Ainsworth M., Levesley J., Marietta M., Light W.) -Oxford: Clarendon Press, 1997.

64. Ben Belgacem F., Maday Y. The mortar element method for three dimensional finite elements. // Modélisation Mathématique et Analyse Numérique (M2AN), 1997, V.31, p.289-302.

65. Ben Belgacem F. The mortar finite element method with Lagrange multipliers. // Numer. Math., 1999, V.84, p.173-197.

66. Benzi M., Golub G., Liesen J. Numerical Solution of Saddle Point Problems. // Acta Numerica, 2005, V.14, p.1-137.

67. Bernardi C, Maday Y., Patera A. A new nonconforming approach to domain decomposition: The mortar element method. // Nonlinear partial differential equations and their applications. ( Brezis H. et al., eds.). -Paris, 1994, p.13-51.

68. Bell J., Trangenstein J., Shubin G. Conservation Laws of Mixed Type Describing Three-Phase Flow in Porous Media. // SIAM J.Appl.Math., V.46, 1986, p.1000-1017.

69. Behie G., Vinsome P. Block iterative methods for fully implicit reservoir simulation. // Soc. of Pet. Eng. J., 1982, V.22, p.658-668.

70. Bernardi C., Maday Y. Mesh adaptivity in finite elements by the mortar method. // Revue européenne des éléments finis, 2000, V.9, p.451-465.

71. Blum H., Lisky S., Rannacher R. A Domain Splitting Algorithm for Parabolic Problems. // Computing, 1992, V.49, p.11-23.

72. Boglaev I. A monotone Schwarz algorithm for a semilinear convection-diffusion problem. // J.Numer.Math., 2004, V.12, p.169-191.

73. Boglaev I. Monotone Schwarz iterates for a semilinear parabolic convection-diffusion problem. // J. of Computational and Applied Mathematics, 2005, V.183, p.191-209.

74. Borouchaki H., George P.-L., Hecht F., Laug P., Mohammadi B., Saltel E. Mailleur bidimensionnel de Delaunay gouverne par une carte de metriques. Rapp. Rech.2760. Inst. Nat. Rech. Informat. Automat. -Paris: INRIA, 1995.

75. Boursier I., Tromeur-Dervout D., Vassilevski Yu. Parallel solution of Mixed Finite Element / Spectral Element systems for convection-diffusion equations on non-matching grids.// Preprint CDCSP-03-01, -Lyon: Université Claude Bernard Lyon 1, 2003.

76. Braess D., Deuflhard P., Lipnikov K. A subspace cascadic multigrid method for mortar elements.// Computing, 2002, V.69, No.3, p.205-225.

77. Bramble J., Pasciak J., Schatz A. The construction of preconditioners for elliptic problems by substructuring, IV. // Math.Comp., 1989, V.53, p.1-24.

78. Bramble J., Pasciak J., Vassilev A. Analysis of non-overlapping domain decomposition algorithms with inexact solves. // Math. Comp., 1998, V.67, p.1-19.

79. Bramble J., Pasciak J., Vassilev A. Analysis of the inexact Uzawa algorithm for saddle point problems. // SIAM J.Numer.Anal., 1997, V.33, No.4, p.1072-1092.

80. Bramble J., Pasciak J. A preconditioning technique for indefinite systems resulting from mixed approximations of elliptic problems. // Math.Comp., 1988, V.50, No.181, p.1-17.

81. Bramble J., Pasciak J., Xu J. Parallel multilevel preconditioners. // Math.Comp., 1990, V.55, p.1-22.

82. Brezina M. Robust iterative solvers on unstructured meshes. PhD thesis, -Denver: University of Colorado at Denver, 1997.

83. Brezina M., Cleary A., Falgout R., Henson V., Jones J., Manteuffel T., McCormick S., Ruge J. Algebraic multigrid based on element interpolation (AMGe). // SIAM J. Sci. Comp., 2000, V.22, p.1570-1592.

84. Brezzi F., Fortin M. Mixed and Hybrid Finite Element Methods. -N.Y.: Springer-Verlag, 1991.

85. Brown P., Saad Y. Hybrid Krylov methods for nonlinear systems of equations. // SIAM J. Sci.Statist.Cornput., 1990, V.ll, p.450-481.

86. Buscaglia G., E. Dari E. Anisotropic mesh optimization and its application in adaptivity. // Inter. J. Numer. Meth. Engng., 1997, V.40, p.4119-4136.

87. Buscaglia G., Agouzal A., Ramirez P., Dari E. On Hessian recovery and anisotropic adaptivity. // Proceedings of ECCOMAS 98, -Athens. 1998, p.403-407.

88. Chan T., Smith B., Zou J. Overlapping Schwarz methods on unstructured meshes using non-matching coarse grids. // Numer. Math., 1996, V.73, No.2, p.149-167.

89. Buzbee B., Golub G., Nielson C. On direct methods for solving Poisson's equations. // SIAM J.Numer.Anal., 1970, V.7, p.627-656.

90. Carey G.F. Computational grids: generation, adaptation, and solution strategies. // Washington: Taylor & Francis, 1997.

91. Chan T., Wan W. Analysis of projection methods for solving linear systems with multiple right-hand sides. // SIAM J. Sci. Comput., 1997, V.18, No.6, p.1698-1721.

92. Chartier T., Falgout R., Henson V., Jones J., Manteuffel T., McCormick S., Ruge J., Vassilevski P. Spectral AMGe. // SIAM J. Sci. Comp., 2004, V.25, No.l, p.1-26.

93. Chen Z., Douglas J. Prismatic mixed finite elements for second order elliptic problems. // Calcolo, 1989, V.26, p.135-148.

94. Chen Z. Equivalence between multigrid algorithms for nonconforming and mixed methods for second-order elliptic problems. // East-West J.Numer.Math., 1996, V.4, p.1-33.

95. Chen Z. On the convergence of Galerkin-multigrid methods for nonconforming finite elements. // East-West J.Numer.Math., 1999, V.7, p.79-104.

96. Chorin A. Numerical solution of the Navier-Stokes equations.// Math. Comp. 1968, V.22, p.745-762.

97. Dendy J. Multigrid methods for petroleum reservoir simulation on SIMD machines. 11 SPE 25243, Proceedings of SPE 1993 Resrv.Simul.Symp., New Orleans, LA, February 28-March 3, 1993, p.97-104.

98. Castro Diaz M., Hecht F., Mohammadi B. New progress in anisotropic grid adaptation for inviscid and viscous flows simulations // Rapp. Rech.2671. Inst. Nat. Rech. Informat. Automat. -Paris: INRIA, 1995.

99. Dompierre J., Vallet M.-G., Fortin M., Habashi W., Ait-Ali-Yahia D., Tam A., Edge-based mesh adaptation for CFD. Technical Report R95-73, CERCA, 1995.

100. Dryja M., 0. Widlund 0. Domain decomposition algorithms with small overlap. // SIAM J. Sci.Comput., 1994, V.15, No.3, p.604-620.

101. Dryja D., Sarkis M., Widlund 0. Multilevel Schwarz methods for elliptic problems with discontinuous coefficients in three dimensions. // Numer.Math., 1996, V.72, p.313-348.

102. Dryja M., Smith B., Widlund 0. Schwarz analysis of iterative substructuring algorithms for elliptic problems in three dimensions. // SIAM J. Numer.Anal., 1994, V.31, p.1662-1694.

103. Dyadechko V., Iliash Yu., Vassilevski Yu. Structuring preconditioners for unstructured meshes. // Russ.J. Numer.Anal. Math.Modelling, 1996, V.ll, p.139-154.

104. Dyadechko V., Lipnikov K., Vassilevski Yu. Hessian based anisotropic mesh adaptation in domains with discrete boundaries. // Russ.J.Numer.Anal.Math.Modelling, 2005, V.20, No.4, p.391-402.

105. Edwards H. A parallel multilevel-preconditioned GMRES solver for multi-phase flow models in the Implicit Parallel Accurate Reservoir Simulator (IPARS). // TIC AM Report 98-04, The University of Texas at Austin. -Austin TX: TICAM, 1998.

106. Eisenstat S., Walker H. Globally convergent inexact Newton methods. // SIAM J. Optim., 1994, V.4, p.393-422.

107. Eisenstat S., Walker H. Choosing the forcing terms in an inexact Newton method. // SIAM J. Sci. Comput., 1996, V.17, p.16-32.

108. Erhel J., Guyomarc'h F. An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems. // SIAM J. Matrix Anal. Appl., 2000, V.21, No.4, p.1279-1299.

109. Farhat C. A simple and efficient automatic FEM domain decomposer. // Computers &; Structures, 1988, V.28, No.5, p.579-602.

110. Farhat C., Roux F. A method of finite element tearing and interconnecting and its parallel solution algorithm. // Int. J. Numer. Meths. Engng., 1991, V.32, p.1205-1227.

111. Feistauer M. Mathematical Methods in Fluid Dynamics. -Harlow, England: Longman Scientific &; Technical, 1993.

112. Fischer P. Projection techniques for iterative solution of Ax = b with successive right-hand sides. // Comput. Methods Appl. Mech. Engrg., 1998, V.163, No.1-4, p.193-204.

113. Fortin M., Vallet M.-G., Dompierre J., Bourgault Y., Habashi W.G. Anisotropic mesh adaptation: theory, validation and applications // Computational Fluid dynamics. John Wiley k Sons Ltd., 1996. p.174-180.

114. Foster I., Designing and building parallel programs. -N.Y.: Addison-Wesley Publishing Company, 1995.

115. Garbey M., Kuznetsov Yu. Parallel Schwarz Algorithm for Equations with Singular-Perturbed Convection-Diffusion Operator. Preprint 238, C.N.R.S.U.M.R. 5585, Université Claude Bernard Lyon 1. -Lyon: UCBL, 1996.

116. Garbey M., Kuznetsov Yu., Vassilevski Yu. Parallel Schwarz method for a convection-diffusion problem.// SIAM J.Sci.Comp., 2000, V.22, No.3, p.891-916.

117. Garbey M., Kuznetsov Yu., Vassilevski Yu. Parallel Schwarz algorithm for a problem with a singular-perturbed convection-diffusion operator. // Preprint CDCSP-97-08, CDCSP, Université Claude Bernard Lyon 1. -Lyon: UCBL, 1997

118. Garbey M., Vassilevski Yu. A parallel solver for unsteady incompressible 3D Navier-Stokes equations.// Parallel Computing, 2001, V.27, No.4, p.363-389.

119. Garbey M., Kuznetsov Yu., Vassilevski Yu. On parallel solution of singularly perturbed convection-diffusion problems. // Proceedings of the Third ECCOMAS Conference on Numerical Methods in Engineering. John Wiley k Sons, 1998, V.2, p.245-250.

120. Garimella R., Swartz B. Curvature estimation for unstructured triangulations of surfaces. Los Alamos Report LA-UR-03-8240, http://math.lanl.gov/Research/Publications/Docs/garimella-2003-curvature.pdf.

121. George P.L., Borouchaki H. Delaunay Triangulation and Meshing. -Paris: Editions Hermes, 1998.

122. Gerardo-Giorda L., Le Tallec P., Nataf F. A Robin-Robin preconditioner for strongly heterogeneous advection-diffusion problems. // Comput. Methods Appl. Mech. Engrg., 2004, V.193, No.9-11, p.745-764.

123. Girault V., Raviart P. Finite elements methods for Navier-Stokes equations. Series in Computational mathematics, 5. -N.Y.: Springer, 1986.

124. Globisch G., Nepomnyaschikh S. The hierarchical preconditioning on unstructured grids. // Computing, 1998, V.61, p.307-330.

125. Glowinski R., Pironneau O., Numerical methods for the first biharmonic equation and for the two-dimensional Stokes problem. // SIAM Rev., 1979, V.21, p.167-212.

126. Gopalakrishnan J., Pasciak J. Multigrid for the Mortar Finite Element Method. // SIAM J. Numer. Anal., 2000, V.37, No.3, p.1029-1052.

127. Gosselet P., Rey Ch. On a selective reuse of Krylov subspaces in Newton-Krylov approaches for nonlinear elasticity. // Domain decomposition methods in science and engineering. -Mexico: Natl. Auton. Univ. Mex., 2003, p.419-426.

128. Graham I., Hagger M. Unstructured additive Schwarz-CG method for elliptic problems with highly discontinuous coefficients. // SIAM J.Sci.Comp., 1999, V.20, p.2041-2066.

129. Grisvard P. Boundary Value Problems in Non-Smooth Domains. -London: Pitman, 1985.

130. Gropp W., Curfman Mclnnes L., B.Smith, Using the scalable nonlinear equations solvers package. Tech.report ANL/MCS-TM-193, -Argonne, IL: Argonne National Laboratory, 1995.

131. Golub G., Van Loan C. Matrix Computations. -North Oxford Academic Publishers Ltd, 1986.

132. Guermond J. Some implementations of projection methods for Navier-Stokes equations. // Model. Math. Anal. Numer. (M2AN), 1996, V.30, p.637-667.

133. Gupta A. Improved symbolic and numerical factorization for unsymmetric sparse matrices. // SIAM J.Matr.Anal.Appl., 2002, V.24, p.529-552.

134. Gunzburger M. Finite element methods for viscous incompressible flows: a guide to theory, practice, and algorithms. -San Diego-London: Academic Press, 1989.

135. Hackbusch W. Iterative solution of large sparse systems of equations. -N.Y.: SpringerVerlag, 1994.

136. Hackbusch W. Multigrid methods and applications. -Berlin: Springer, 1985.

137. Hakopian Yu., Kuznetsov Yu. Algebraic multigrid/substructuring preconditioners. // Sov. J. Numer. Anal. Math. Modelling, 1991, V.6, p.453-483.

138. Hoppe R., Iliash Yu., Kuznetsov Yu., Vassilevski Yu., Wohlmuth B. Analysis and parallel implementation of adaptive mortar element methods. // East-West J. Numer. Math., 1998, V.6, No.3, p.223-248.

139. Hoppe R., Petrova S., Vassilevski Yu. Adaptive grid refinement for computation of the homogenized elasticity tensor. // Proceedings of 4th Intern. Conf. on Large-Scale Sei. Comput. LSSC-03, Lect.Notes Comput.Sci., Springer, 2004, V.2907, p.371-378.

140. Iliash Yu., Kuznetsov Yu., Vassilevski Yu. Efficient parallel solvers for two dimensional potential flow and convection-diffusion problems on nonmatching grids. // Report No. 357, Universität Augsburg. -Augsburg, Germany, 1996.

141. Jenkins E., Kees C., Kelley C., Miller C. An aggregation-based domain decomposition preconditioner for groundwater flow. // SIAM J. Sci. Comp., 2001, V.25, p.430-441.

142. Johnson C., Schatz A., Wahlbin L. Crosswind smear and pointwise errors in streamline diffusion finite element methods. // Math.Comp., 1987, V.49, p.25-38.

143. Keshtiban I., Belblidia F., Webster M. Compressible flow solvers for Low Mach number flows a review. // Technical report CSR2, Institute of Non-Newtonian Fluid Mechanics, University of Wales. -Swansea, UK, 2004.

144. Kincaid D., Respess J., Young D., Grimes R. Algorithm 586 ITPACK 2C: A fortran package for solving large sparse linear systems by adaptive accelerated iterative methods. // ACM Transactions on Mathematical Software, 1982, V.8, p.302-322.

145. Kim C., Lazarov R., Pasciak J., Vassilevski P. Multiplier spaces for the mortar finite element method in three dimensions. // SIAM J. Numer. Analysis, 2001, V.39, p.519-538.

146. Kim J., Moin P. Application of a fractional step method to incompressible Navier-Stokes equations. // J. Comput. Phys., 1985, V.59, p.308-323.

147. Kiss B., Krebsz A. A sparse representation of the Hnorm in finite element spaces. Preprint N 1, Dep. of Math., College Gyor, Hungary, 1995.

148. Kaporin I. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + ^/-decomposition. // Numer.Linear Algebra Appl., 1998, V.5, p.483-509.

149. Klie H. Krylov-secant methods for solving large scale systems of coupled nonlinear parabolic equations. PhD Thesis, Rice University, Houston, 1996.

150. Kunert G. Aposteriori error estimation for anisotropic tetrahedral and triangular finite element meshes. Dissertation, 1999, TU Chemnitz.

151. Kuznetsov Yu. Matrix computational processes in subspaces. 11 Comp.Math. Appl.Sci.Eng. -Amsterdam: North-Holland, 1984, V.6, p.15-31.

152. Kuznetsov Yu. New Algorithms for Approximate Realization of Implicit Difference Schemes. // Sov. J. Num. Anal. Math. Modelling, 1988, V.3, p.99-114.

153. Kuznetsov Yu. Overlapping domain decomposition methods for finite element problems with singular perturbed operators. // DDM for PDE's, Proc. of 4th Int.Symp. (Glowinski R. et al., eds), -Philadelphia: SIAM, 1991, p.223-242.

154. Kuznetsov Yu. Efficient iterative solvers for elliptic problems on nonmatching grids. // Russ. J. Numer. Anal. Math. Modelling, 1995, V.10, No.3, p.187-211.

155. Kuznetsov Yu., Wheeler M. Optimal order substructuring preconditioners for mixed finite element methods on nonmatching grids. // East-West J. Numer. Math., 1995, V.3, p.127-143.

156. Kuznetsov Yu., Rossi T. Fast direct method for solving algebraic systems with separable symmetric band matrices. // East-West J. Numer.Math., 1996, V.4, p.53-68.

157. Kuznetsov Yu., Vassilevski Yu. An interface preconditioner for the mortar element method. // Numerical Mathematics and Advanced applications, Proc. 3d European Conference (Neitaanmaki P. et al., eds.), World Scientific, 2000, p.753-761.

158. Kuznetsov Yu., Manninen P., Vassilevski Yu. On numerical experiments with Neumann-Neumann and Neumann-Dirichlet domain decomposition preconditioners. Technical report, University of Jyvaskyla, Finalnd, 1993.

159. Lacroix S., Vassilevski Yu., Wheeler M., Wheeler J., Iterative solution methods for modeling multiphase flow in porous media fully implicitly.// SIAM J.Sci.Comp., 2003, V.25, No.3, p.905-926.

160. Lacroix S., Vassilevski Yu., Wheeler M. Decoupling preconditioners in the Implicit Parallel Accurate Reservoir Simulator (IPARS).// Numer. Linear Algebra Appl., 2001, V.8, No.8, p.537-549.

161. Lehoucq R., Sorensen D., Yang C. ARPACK Users' guide: solution of large-scale eigenvalue problems with implicitly restarted arnoldi methods. //Software, Environments, and Tools V.6. -SIAM, 1998.

162. Lipnikov K., Vassilevski Yu. Parallel adaptive solution of 3D boundary value problems by Hessian recovery.// Comp.Methods Appl.Mech.Engnr., 2003, V1.192, No.11-12, p.1495-1513.

163. Lipnikov K., Vassilevski Yu. Parallel adaptive solution of the Stokes and Oseen problems on unstructured 3D meshes. // Proceedings of Int.Conf. Parallel CFD 2003, Elsevier, 2004, p.153-162.

164. Lipnikov K., Vassilevski Yu. On control of adaptation in parallel mesh generation. // Engrg. with Comput., 2004, V.20, No.3, p.193-201.

165. Lipnikov K., Vassilevski Yu. Error estimates for Hessian-based mesh adaptation algorithms with control of adaptivity. // Proc. of 13th Int. Meshing Roundtable, Williamsburg, Virginia, 2004, p.345-351

166. Lipnikov K., Vassilevski Yu. On a parallel algorithm for controlled Hessian-based mesh adaptation.// Proc. of 3d Conf. Appl.Geometry, Mesh Generation and High Performance Computing, -M: Comp.Center RAS, 2004, V.l, p. 155-166.

167. Lipton R., Rose D,, Tarjan R. Generalized nested dissection. // SIAM J.Numer.Anal., 1979, V.16, p.346-358.

168. Lisejkin V. Grid generation methods. -Berlin: Springer, 1999.

169. Mandel J. Hybrid domain decomposition with unstructured subdomains. // Domain Decomposition Methods in Science and Engineering, Proc. of the 6th Int. Conf. on DDM, Contemporary Mathematics, AMS, 1994, V.157, p.103-112.

170. Mandel J. Balancing domain decomposition. // Comm. Numer. Meth. Engrg., 1993, V.9, p.233-241.

171. Mandel J., Bresina M. Balancing domain decomposition for problems with large jumps in coefficients. // Math.Comp., 1996, V.65, p.1387-1401.

172. Marchuk G., Kuznetsov Yu., Matsokin A. Fictitious domain and domain decomposition methods. // Sov.J.Numer.Anal. and Math.Modelling, 1986, V.l, No.l, p.3-35.

173. Marion M., Temam R. Navier-Stokes equations: theory and approximation. // Handbook of numerical analysis, V.VI. -Amsterdam: Elsevier, 1998, p.503-689.

174. Matsokin A., Nepomnyaschikh S. A Schwarz alternating method in a subspace. // Soviet Mathematics, 1985, V.29, No.10, p.78-84.

175. Mclvor A., Valkenburg R. A comparison of local surface geometry estimation methods. // Machine Vision and Applications, 1997, V.10, p.17-26.

176. Meijerink J., Van der Vorst H. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. // Math.Comp., 1977, V.31, p.148-162.

177. Meyer M., Lee H., Barr A., Desbrun M. Generalized barycentric coordinates on irregular polygons. // Journal of Graphics Tools, 2002, V.7, No.l, p.13-22.

178. Mokhtarian F., Khalili N., Yuen P. Curvature computation of free-form 3-D meshes at multiple scales. // Computer Vision and Image Understanding, 2001, V.83, p.118-139.

179. Nepomnyaschikh S. Mesh theorems of traces, normalizations of function traces and their inversion.// Sov.J.Numer.Anal.Math.Modelling, 1991, V.6, p.223-242.

180. Nepomnyaschikh S. Fictitious space method on unstructured meshes. // East-West J.Numer.Math., 1995, V.3, p.71-79.

181. Nepomnyaschikh S. Finite element trace theorems for parameter dependent sobolev spaces. // Numerical Mathematics and Advanced applications, Proc. of 3d European Conference (Neitaanmaki P. et al., eds.), World Scientific, 2000, p. 31-41.

182. Nepomnyaschikh S. Application of domain decomposition to elliptic problems on with discontinuous coefficients. // Proc. of the 4th Int. Symp. on DDM for PDE's (Glowinski R. et al., eds.). -Philadelphia, PA: SIAM, 1991, p.242-251.

183. Oliphant T. An implicit numerical method for solving two-dimensional time-dependent diffusion problems.// Quart.Appl.Math., 1961, V.19, p.221-229.

184. O'Neil J., Szyld D. A block ordering method for sparse matrices. // SIAM J. Sci. Stat. Comput., 1990, V.ll, No.5, p.811-823.

185. Niijima K. Pointwise error estimates for a streamline diffusion finite element scheme. // Numer.Math., 1990, V.56, p.707-719.

186. Olshanskii M., Reusken A. Convergence analysis of a multigrid solver for a finite element method applied to convection-dominated model problem. // SIAM J.Num.Anal., 2004, V.43, p.1261-1291.

187. Peraire J., Morgan K., Peiro J. Unstructured mesh methods for CFD // I.C. Aero Report 90-04. Imperial College. UK. 1990.

188. Pernice M., Walker H. NITSOL: a Newton iterative solver for nonlinear systems. // SIAM J. Sci. Comput., 1998, V.19, p.302-318.

189. Peskin C., The immersed boundary method.// Acta Numerica, 2002, p.1-39.

190. Pironneau O., Finite element methods for fluids. -Chichester: J. Wiley & Sons, 1989.

191. Pothen A., Simon H., Liou K. Partitioning sparse matrices with eigenvectors of graphs. // SIAM J. Math. Anal. Appl., 1990, V.ll, p.430-452.

192. Proceedings of the 15th Symposium on Reservoir Simulation, Society of Petroleum Engineers Inc., Richardson TX, 1999.

193. Prudhomme S., Oden J., On goal-oriented error estimation for elliptic problems: Application to the control of pointwise errors. // Comput. Methods Appl. Mechan. Engrg., 1999, V.176, p.313-331.

194. Quarteroni A., Valli A. Domain decomposition methods for partial differential equations. Oxford Science Publications, 1999.

195. Rixen D., Farhat C., Tezaur R., Mandel J. Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver. // Int. J. Numer. Meth. Engrg., 1999, V.46, p.501-534.

196. Rannacher R., Zhou G. Analysis of domain splitting method for non-stationary convection-diffusion problems. // East-West J.Numer.Math., 1994, V.2, p.151-172.

197. Rathinam M., Petzold L. Dynamic iteration using reduced order models: a method for simulation of large scale modular systems. // SIAM J. Numer.Anal., 2002, V.40, p.1446-1474.

198. Rathinam M., Petzold L. A new look at proper orthogonal decomposition. // SIAM J. Numer.Anal., 2003, V.41, p.1893-1925.

199. Raviart R, Thomas J. A mixed finite element method for 2-nd order elliptic problems. 11 Mathematical Aspects of Finite Element Methods (Galligani I. and Magenes E., eds.). Lecture Notes in Mathematics, V.606. -N.Y.:Springer-Verlag, 1977, p.192-315.

200. Rippa S. Long and thin triangles can be good for linear interpolation // SIAM J. Numer. Anal, 1992, V.29, p.257-270.

201. Risler F., Rey Ch. Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale nonlinear problems. // Numer. Algorithms, 2000, V.23, No.l, p.1-30.

202. Rivara M. Mesh refinement processes based on the generalized bisection of simplexes. // SIAM J.Numer.Anal., 1984, V.21, p.604-613.

203. Rivara M. Selective refinement/derefinement algorithms for sequences of nested triangulations. // Int.J.Numer.Meth.Engrg., 1989, V.28, p.2889-2906.

204. De Roeck Y.-H., Le Tallec P. Analysis and test of a local domain decomposition preconditioner. // Proc. of the 4th Int. Symp. on DDM for PDE's (Glowinski R. et al., eds.). -Philadelphia, PA: SIAM, 1991, p.112-128.

205. Rossi T., Toivanen J. A parallel fast direct solver for block tridiagonal systems with separable matrices of arbitrary dimension. // SIAM J.Sci.Comp., 1999, V.20, p.1778-1793.

206. Rossi T. Fictitious domain methods with separable preconditioners. PhD thesis, Dept. of Mathematics, University of Jyvaskyla, Finland, 1995.

207. Rosar M., Peskin C. Fluid flow in collapsible ellastic tubes: a three-dimensional numerical model. // New York J. Math., 2001, V.7, p.281-302.

208. Saad Y. Iterative methods for sparse linear systems, Second Edition. -Philadelphia, PA: SIAM, 2003.

209. Saad Y. A flexible inner-outer preconditioned GMRES algorithm. // SIAM J.Sci.Comput., 1993, V.14, p.461-469.

210. Schefer M., S.Turek S. Benchmark computations of laminar flow around a cylinder. // Flow simulation with high-performance computers II (Hirschel E., eds.), Notes on numerical fluid mechanics, V.52, 1996. -Braunschweig: Vieweg, p.547-566.

211. Siebert K. An a posteriori error estimator for anisotropic refinement. // Numer. Math., 1996, V.73, p.373-398.

212. Smith B., Bjorstad P., Gropp W. Domain decomposition: parallel multilevel methods for elliptic partial differential equations. // Cambridge University Press, 1996.

213. Stefanica D., Klawonn A. The FETI method for mortar finite elements. // Proc. of 11th Int. Conf. on DDM (Lai C.-H. et al.,eds.), DDM.org, 1999, p.121.129.

214. Stiiben K. A review of algebraic multigrid. // J.Comp. and Appl.Mathematics, 2001,1. V.128, No.1-2, p.281-309.

215. Stiiben K. Algebraic multigrid (AMG): experiences and comparisons. // Appl. Math. Comput., 1983, V.13, p.419-452.

216. Temam R. Sur l'approximation des équations de Navier-Stokes. C.R.Acad.Sci. Paris, Série A, 1966, V.262, 219-221.

217. Tomasset F. Implementation of finite element for Navier-Stokes equations. // N.Y.: Springer-Verlag, 1981.

218. Trangenstein J., Bell J. Mathematical structure of the black-oil model for petroleum reservoir simulation. // SIAM J.Appl.Math., 1989, V.49, p.749-783.

219. Tromeur-Dervout D., Vassilevski Yu. Choice of initial guess in iterative solution of series of systems.// Submitted to J.Comp.Phys.

220. Tromeur-Dervout D., Vassilevski Yu. POD acceleration of fully implicit solvers for unsteady nonlinear problems and its GRID applications. // Advances in Engineering Software, 2005, to appear.

221. Tromeur-Dervout D. Résolution des Equations de Navier-Stokes en Formulation Vitesse Tourbillon sur Systèmes multiprocesseurs à Mémoire Distribuée. PhD Thesis, Univ. Paris1. VI, 1993.

222. Turek S. Efficient solvers for incompressible flow problems: an algorithm approach in view of computational aspects. -Berlin-Heidelberg: Springer-Verlag, 1999.

223. Tyrtyshnikov E. Mosaic ranks and skeletons. // Lecture Notes in Comput. Sci., V.1196. -Berlin: Springer, p.505-516.

224. Tyrtyshnikov E. Mosaic-skeleton approximations. // Calcolo, 1998, V.33, p.47-57.

225. Tyrtyshnikov E. Incomplete cross approximation in the mosaic-skeleton method. // Computing, 2000, V.64, p.367-380.

226. Tyrtyshnikov E., Vassilevski Yu. A mosaic preconditioner for a dual schur complement.// Numerical Mathematics and Advanced Applications, Proc. of ENUMATH'01. -Milano: Springer-Verlag Italia, 2003, p.867-880.

227. Vangk P., Brezina M., Mandel J. Convergence of algebraic multigrid based on smoothed aggregation. // Numer.Math., 2001, V.88, p.559-579.

228. Varga R. Factorizations and normalized iterative methods. // Boundary problems in differential equations. -Madisson WI: Univ. of Wisconsin Press, 1960, p.121-142.

229. Vassilevski Yu. A hybrid domain decomposition method based on aggregation.// Numer. Linear Algebra Appl., 2004, V.ll, p.327-341.

230. Vassilevski Yu. A parallel interface preconditioner for the mortar element method in case of jumping coefficients.// Domain Decomposition Methods in Sciences and Engineering, DDM.org, 2001, p.231-240.

231. Vassilevski Yu. Domain decomposition methods and averaging operators for the case of multidomain splitting. // Russ.J.Numer.Anal.Math.Modelling, 1995, V.10, No.2, p.141-148.

232. Vassilevski Yu. A parallel CG solver based on domain decomposition and non-smooth aggregation. // Conjugate gradient algorithms and finite element methods (Proceedings of Int.Conf. 50 years of CG), -Berlin-Heidelberg: Springer-Verlag, 2004, p.93-102.

233. Vassilevski P. Fast algorithm for solving a linear algebraic problem with separable variables. // C.r.Acad.Bulg.Sci., 1984, V.37, p.305-308.

234. Verfurth R. A review of a posteriori error estimation and adaptive mesh-refinement techniques. -Stuttgart: Wiley-Teubner, 1996.

235. Wallis J., Kendall R., Little T. Constrained Residual Acceleration of Conjugate Residual Methods. // SPE 13536, Proc. of the SPE 1985 Resrv.Simul.Symp., (Dallas, Texas), 1985, p.415-428.

236. Watts J. An iterative matrix solution method suitable for anisotropic problems.// Soc.Pet.Eng.J., 1971, p.47-51.

237. Watts J. A method of improving line successive overrelaxation in anisotropic problems a theoretical analysis. // Soc.Pet.Eng.J., 1973, p.105-118.

238. Watts J. A conjugate gradient truncated direct method for the iterative solution of the reservoir simulation pressure equation. // Soc.Pet.Eng.J., 1981, V.21, p.345-353.

239. Watts J. A compositional formulation of the pressure and saturation equations. // SPE Reservoir Engineering, 1986, V.l, p.243-252.

240. Wheeler J., Smith R. Reservoir simulation on a hypercube. // SPE 19804, Proc. of the SPE 1989 Ann.Tech.Conf., (San Antonio, Texas), 1989.

241. Wheeler J. IPARS user's manual. Report CSM, TICAM, The University of Texas at Austin, 2000.

242. Wheeler M., Yotov I. Multigrid on the interface for mortar mixed finite element methods for elliptic problems. // Comp.Meth.Appl.Mech.Engng., 2000, V.184, p.287-302.

243. Williams R. Performance of dynamic load balancing algorithms for unstructured mesh calculations. // Concurrency: Practice and experience, 1992, V.3, p.457-481.

244. Wohlmuth B. A residual based error estimator for mortar finite element discretizations. // Numer. Math., 1999, V.84, p.143-171.

245. Wohlmuth B. A multigrid method for saddle point problems arizing from mortar finite element discretizations. // Electronic Transactions on Numerical Analysis, 2000, V.ll, p.43-54.

246. Young D. Iterative solution of large linear systems. -N.Y.: Academic Press, 1971.

247. Zavattieri P., Dari E., Buscaglia G. Optimization strategies in unstructured mesh generation. // Inter. J. Numer. Meth. Engng., 1996, V.39, p.2055-2071.

248. Handbook of grid generation (Tompson J., Soni B., Weatherill N., eds.). -Boca Raton, FL: CRC Press, 1999.275. http://www.ddm.org/Bibtex/bibtex.html