Содержание к диссертации
Введение
Глава 1 Блочные параллельные методы для решения неструктурированных систем 34
1.1 Последовательные технологии решения неструктурированных систем . 35
1.1.1 Методы факторизации 35
1.1.2 Методы неполной факторизации 38
1.1.3 Алгебраические многосеточные методы 40
1.1.4 Структурированные переобуславливатели для неструктурированных сеток 41
1.2 Решение неструктурированных систем с матрицами тензорного ранга 2 . 43
1.3 Блочные переобуславливатели и переупорядочивание для задач с анизотропными коэффициентами 46
1.4 Параллельный метод агрегирования и его применения 52
1.4.1 Метод агрегирования 53
1.4.2 Параллелизация метода агрегирования 61
1.4.3 Диффузионные задачи с гетерогенными коэффициентами и/или большим числом подобластей 63
1.4.4 Параллелизация на адаптивных неструктурированных сетках . 65
Выводы главы 1 68
Глава 2 Параллельные адаптивные технологии на симплициальных анизотропных сетках 69
2.1 Оптимальные и квази-оптимальные сетки и их свойства 70
2.1.1 Оптимальные сетки 70
2.1.2 Квази-оптимальные сетки 78
2.2 Последовательный адаптивный алгоритм построения квази-оптимальных сеток 84
2.2.1 Восполнение сеточного гессиана 84
2.2.2 Генерация сеток, квази-равномерных в заданной метрике 92
2.2.3 Адаптивный алгоритм 95
2.3 Управление адаптацией 98
2.4 Особенности адаптации для областей с дискретной границей 106
2.5 Параллельный адаптивный алгоритм для приближенного решения краевых задач 113
2.5.1 Параллельная генерация сеток, квази-равномерных в заданной метрике 113
2.5.2 Параллельное адаптивное решение трехмерных краевых задач . 116
2.5.3 Влияние управления адаптацией 124
Выводы главы 2 128
Глава 3 Блочные параллельные методы решения неконформных конечно-элементных систем 129
3.1 Макро-гибридная постановка на нестыкующихся сетках 130
3.2 Интерфейсные переобуславливатели 137
3.2.1 Переобуславливатели с внутренним итерационным процессом для множителей Лагранжа 138
3.2.2 Дирихле-Дирихле переобуславливатель 143
3.2.3 Численные эксперименты с параллельными решателями 146
3.3 Параллельные адаптивные технологии на нестыкующихся сетках 151
3.4 Мозаично-скелстонные переобуславливатели для множителей Лагранжа . 156
3.5 Параллельный многосеточный метод для неконформных конечных элементов 165 Выводы главы 3 172
Глава 4 Блочные параллельные технологии на трехмерных прямоугольных сетках 173
4.1 Параллельный двухуровневый метод Шварца для сингулярно возмущенного уравнения конвекции-диффузии 174
4.1.1 Двухуровневый метод Шварца для постоянного вектора переноса . 175
4.1.2 Двухуровневый метод Шварца для переменного вектора переноса . 182
4.1.3 Адаптивный алгоритм разбиения для случая переменного вектора переноса 192
4.2 Некоторые приложения параллельного двухуровневого метода Шварца . 193
4.2.1 Трехмерное уравнение конвекции-диффузии 193
4.2.2 Трехмерное уравнение конвекции-диффузии-реакции 196
4.2.3 Проекционный алгоритм для нестационарных уравнений Навье-Стокса196
4.2.4 Проекционный алгоритм для уравнения тепло-массопереноса в простейшем химическом реакторе 201
4.3 Адаптивные итерационные технологии для последовательности систем . 204
4.3.1 Последовательность систем с одной матрицей и разными правыми частями 204
4.3.2 Последовательность систем с разными матрицами и правыми частями207
4.3.3 Последовательность нелинейных систем 210
4.4 Блочные параллельные методы для систем уравнений многофазной фильтрации 219
4.4.1 Уравнения многофазной фильтрации, их дискретизация и линеаризация 219
4.4.2 Многосеточный метод с ускоренным разгрублением 221
4.4.3 Алгебраическое выделение уравнения для давления 225
Выводы главы 4 230
Заключение 232
- Блочные переобуславливатели и переупорядочивание для задач с анизотропными коэффициентами
- Параллельная генерация сеток, квази-равномерных в заданной метрике
- Переобуславливатели с внутренним итерационным процессом для множителей Лагранжа
- Проекционный алгоритм для уравнения тепло-массопереноса в простейшем химическом реакторе
Введение к работе
Настоящая диссертация посвящена рассмотрению некоторых современных вычислительных технологий нахождения приближенного решения краевых задач с эллиптическими операторами второго порядка. Под вычислительными технологиями здесь понимается совокупность алгоритмов, структур данных, расчетных методик и программных реализаций для решения вычислительных задач на вычислительных системах. Главной целью вычислительных технологий является получение ответа в нужной форме, к заданному сроку и при известных ограничениях на доступ к машинным и человеческим ресурсам [8] . Основные черты рассматриваемых здесь технологий — параллельность и эффективность.
Параллельность вычислительных технологий вызвана необходимостью решать задачи настолько большой размерности, что даже новейшие однопроцессорные компьютеры не в состоянии обеспечить требуемый расчет. Помимо практических ограничений на быстродействие процессора и размеры компьютерной памяти, немаловажно учитывать и экономический аспект вычислений: очевидна нелинейная зависимость цены компьютера от объема памяти и реального быстродействия процессора. В силу этого параллелизация вычислений, во-первых, позволяет решать гораздо большие задачи, и во-вторых, обеспечивает примерно линейную зависимость цены параллельного компьютера с распределенной памятью от его быстродействия и объема памяти. Именно последнее обстоятельство объясняет широкую распространенность параллельных компьютеров с распределенной памятью и естественный запрос на параллельные технологии, позволяющие использовать такие компьютеры. Распределенная память подразумевает разбиение данных на блоки, обрабатываемые каждым процессором, поэтому блочность алгоритмов характерна для большинства параллельных методов и технологий [18].
Эффективность вычислительных технологий связана с разумным распределением и использованием вычислительных ресурсов, для достижения заданной точности расчетов минимальными вычислительными затратами, или, что равнозначно, для повышения точности расчетов на заданной вычислительной системе. Одним из самых мощных средств повышения эффективности технологии является ее адаптация к конкретной зада-
че. Адаптивность технологии может пониматься в очень широком смысле: от возможности подстраивать выполнение базовых арифметических операций к конкретной архитектуре ЭВМ, до адаптации самого вычислительного алгоритма к особенностям конкретной задачи. В данной работе рассматриваются два вида адаптивности: адаптивное построение расчетной сетки, и адаптация алгоритма решения дискретных задач.
Таким образом, ключевыми инструментами при разработке параллельных и эффективных вычислительных технологий являются блочность и адаптивность, которые проявляются на двух основных этапах решения краевых задач — построении расчетных сеток и дискретизаций, и решении порожденных ими систем.
Общепринятым подходом при разработке современных технологий решений краевых задач математической физики является их замена сеточными уравнениями, обеспечивающая приближенное решение этих задач. Одной из основных составляющих технологии построения систем сеточных уравнений является генерация расчетной сетки. Адаптивно построенные сетки позволяют достигать желаемой точности при приближении исходной краевой задачи на меньшем количестве узлов и элементов сетки. Наиболее широкое освещение вопросов адаптивного приближенного решения краевых задач можно найти в книгах Р.Ферфурта [260], И.Бабушки и Т.Струболиса [59], хотя некоторые ключевые моменты технологии были рассмотрены уже в монографии Л.А.Оганесяна и Л.А.Руховца [33]. Вопросы генерации адаптивных сеток рассматривались в монографиях П.Джорджа и Х.Буручаки [143], В.Д.Лизейкина [193], Г.Кэри [94], а также в сборнике [274]. Системы сеточных уравнений могут решаться как прямыми, так и итерационными методами, однако применение прямых методов ограничено лишь линейными системами с матрицами умеренных порядков или специального вида. Хорошее представление о состоянии теории и практики итерационных методов решения сеточных уравнений можно получить из книг Д.Янга [272], Г.И.Марчука [32], Г.И.Марчука и В.И.Лебедева [30], А.А.Самарского и Е.С.Николаева [37], Е.Г.Дьяконова [22], В.Хакбуша [157], И.Саада [234], Дж.Деммеля [20].
Перейдем к общему обоснованию предлагаемых в диссертационной работе технологических решений.
Использование адаптивных сеток для построения сеточных уравнений зачастую подразумевает решение неструктурированных систем, графы матриц которых не обладают никакой структурой. Аппроксимации на структурированных сетках позволяют использовать эффективные методы решения линейных систем. Так, дискретизации на прямоугольных сетках предполагают лексикографическое упорядочивание узлов и связанных с ними неизвестных в уравнении, что позволяет использовать быстрые прямые методы для решения соответствующих систем [93, 62, 259, 28, 180, 231]. Если сетка представляет
собой подмножество ячеек прямоугольной сетки (как в методе фиктивных компонент), то возможно такое упорядочивание узлов, что быстрые прямые методы применимы в качестве переобуславливателей [197, б]. Использование иерархических локально сгущающихся сеток позволяет использовать многосеточные технологии [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, т.е. иредставимы в виде Лі<8>М2з + МіЛ2з- Здесь Л23, М2з (Лі, Мі) — матрицы жесткости и масс (лампироваииая) на двумерном (соотв.,одномерном) сечении призматической сетки. Для решения таких систем автором предлагается новый метод, обладающий почти линейной арифметической сложностью решения и невысокой сложностью инициализации. Метод представляет собой обобщение быстрого прямого метода [175, 259, 28, 43, 180, 232, 231] на случай неструктурированных матриц Л23 с двумерным графом и использует явный вид матриц А\, Л2з,Мі, М23. Смысл обобщения — замена прямого решения двумерных систем с матрицами Л2з + АМ2з с произвольным параметром А на приближенное итерационное решение с иереобуславливателем, базирующемся на факторизации одной из трех матриц Л23+ЮМ23, Л2з+150М2з, Л2з+1500М23. При этом результирующий переобуславливающий оператор будет нелинейным, и метод из прямого становится быстро сходящимся итерационным (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/minp). Если снабдить блочно-диагональный переобуславливатель коррекцией на каркасном подпространстве, оценка числа обусловленности улучшится то С(р) (1 + у) [92], а в случае гладкого коэффициента р ее можно довести до С (l + у) [118, 238].
Здесь и далее C(z) обозначает константу, зависящую от z, но не зависящую от остальных параметров. Кроме этого, в дальнейшем с, С (с индексами и без индексов) будут обозначать некоторые положительные константы, не зависящие от параметров задачи, а выражение А ~ В будет обозначать спектральную эквивалентность матриц А и В с константами, не зависящими от размеров матриц, или приближенное равенство или пропорциональность величин А и В.
Отметим, что зависимость к от р не является критической, если коэффициент р(х) гладкий в больших подобластях и применяется метод сопряженных градиентов (СГ). Как показано в [150], при сильных скачках число итераций СГ лишь незначительно возрастает на величину, зависящую от геометрической структуры скачков р(х). Однако этот результат не применим для других итерационных методов, а также для гетерогенного коэффициента р(х).
В настоящей работе предлагается метод на основе аддитивного метода Шварца, впервые обеспечивающий теоретическую независимость к от гетерогенного коэффициента диффузии р(х) в зоне перекрывания и оценку к < С (l + у) (Теорема 1.3). В доказательстве теоремы используется гладкость р(х) в подобластях вне зоны перекрывания, хотя это не нужно на практике. Важным практическим достоинством метода является его независимость от гетерогенности р(х) в зоне перекрывания, что позволяет не координировать разрезание на подобласти с линией скачка. Суть метода — специфическое каркасное подпространство, получаемое агрегированием неизвестных внутри подобластей. Агрегирование неизвестных — весьма востребованная технология при разработке числен-
ных методов [86, 164, 253]. Предлагаемый метод, как и метод [164], использует негладкие базисные функции каркасного простраснтва, что отличает его от работ [86, 253]. В отличие от [164], наше каркасное пространство является гораздо большим, поскольку оно включает в себя базисные функции на мелкой сетке, которые не равны нулю в зоне перекрывания. С одной стороны, это делает каркасную (агрегированную) задачу более сложной для решения, с другой стороны, это позволяет исключить зависимость к, от р(х).
Параллелизация решения агрегированной задачи зависит от конкретного приложе
ния. В случае двумерных сеток, уменьшение размерности агрегированной матрицы по
отношению к исходной размерности настолько существенно (N — iV1/2), что позволя
ет использовать прямой метод (факторизации) на каждом процессоре [78, 79]. В случае
трехмерных сеток каркасная задача остается достаточно большой (понижение размерно
сти поэтому ее приближенное решение осуществляется с помощью блочно
го метода верхней последовательной релаксации BSOR или его симметризованного ана
лога. Эффективность подобного решения мотивируется тем, что арифметическая цена
одной итерации BSOR существенно ниже умножения на исходную матрицу жесткости,
а жесткость агрегированной матрицы существенно ниже жесткости исходной матрицы,
что порождает высокую скорость сходимости BSOR. Легкая параллелезуемость метода
является его привлекательной чертой, также как простота его реализации на неструкту
рированных сетках [187,188, 184]. По сути дела, метод является методом "серого ящика",
поскольку помимо элементов исходной матрицы жесткости, требует только разбиение
степеней свободы на непересекающиеся подмножества, обеспечиваемые, например, лю
бым алгоритмом разбиения графа. Отметим, что по формальным признакам предлага
емый алгоритм относится к гибридным методам [194, 238], поскольку его аддитивная
структура на уровне подобластей сочетается с мультипликативной процедурой коррек
ции в каркасном пространстве.
Доступные методы решения неструктурированных систем позволяют использовать полностью неструктурированные симплициальные сетки в адаптивных технологиях приближенного решения краевых задач, которые, в свою очередь, открывают путь к анизотропной адаптации, рассмотренной в главе второй. Как уже отмечалось выше, именно возможность анизотропной адаптации расширяет границы современных адаптивных технологий, поскольку позволяет эффективно аппроксимировать анизотропные решения. В ряде работ [111, ПО, 112, 226], появившихся на рубеже 90-х годов, было показано, что треугольники с тупыми и острыми углами, вытянутые вдоль направления минимальной второй производной некоторой функции, могут быть наилучшими элементами для минимизации ошибки интерполяции. При этом степень вытянутости треугольников и
их малости зависит от матрицы вторых производных (гессиана) функции. На практике управление вытянутостыо элементов задавалось метрикой на основе гессиана функции или его приближения [90, 117, 77, 213, 134]. Предпринимались попытки обобщений на адаптацию к вектор-функциям [116]. Однако первое теоретическое обоснование анизотропной адаптации появилось в работах [11, 53], где были доказаны асимптотические аппроксимационные свойства оптимальных симплициальных сеток и их приближений, квази-оптимальных сеток, которые можно строить на практике.
Технология адаптивного приближенного решения на основе восстановления гессиана сеточного решения породила ряд теоретических вопросов. На некоторые из них даны ответы во второй главе в виде теоретических результатов автора. Так, доказано существование оптимальных сеток при разумных предположениях на непрерывность ошибки интерполирования от координат сеточных узлов и не увеличении этой ошибки при иерархическом измельчении треугольников [12] (Теорема 2.3). Выведены асимптотические оценки ошибки в норме Loo Для оптимальных сеток, в двумерном [11] и трехмерном [53] случаях, предложен их объединенный анализ [10] (Теорема 2.4), а в Теореме 2.7 приведены обобщения (совместно с А.Агузалом) на случай норм Lp, 0 < р < +оо [10]. Рассмотрены свойства приближений оптимальных сеток, которые являются квазиравномерными в метрике, порожденной гессианом функции, и доказана Теорема 2.8 о том, что такие приближения являются квази-оптимальными сетками, ошибка интерполяции на которых асимптотически имеет такую же зависимость от числа симплексов сетки, что и на оптимальных сетках [11, 53].
Восстановление гессиана сеточной функции — важная часть адаптивной технологии. По сути дела, необходимо вычислить матрицу вторых производных у сеточного решения, что на практике осуществляется в рамках конечно-элементных аппроксимаций (за счет перебрасывания одной производной на пробную функцию в формуле Грина). Подобное восстановление не обеспечивает поточечную сходимость восстановленного гессиана к дифференциальному для функций с особенностями, поэтому предположение о поточечной аппроксимации, используемое в излагаемой теории, может не выполняться. С учетом этого автором был предложен новый метод восстановления [54] гессиана у функций с особенностями, который обеспечивает локальную аппроксимацию компонент дифференциального гессиана в L\ норме, что имеет смысл для функций из W?. Здесь и далее W% обозначает соболевское пространство вещественных функций, у которых р-е степени производных порядка не выше q интегрируемы в рассматриваемой области. Отметим, что решения большинства эллиптических краевых задач принадлежат пространству W? в подобластях, где коэффициенты дифференциального оператора — гладкие функции [151]. Это открывает путь к локальной аппроксимации гессиана сеточных функций с
особенностями. Доказанная автором Теорема 2.15 о локальной сходимости восстановленного гессиана к дифференциальному — первый теоретический результат для функций с особенностями и анизотропных сеток.
Генерация квази-равномерных сеток в заданной непрерывной метрике — ключевой элемент рассматриваемой адаптивной технологии [90, 117, 77, 213, 134]. В работе предлагается совместный алгоритм автора и К.Н.Липникова генерации тетраэдральных сеток, являющийся обобщением двумерного алгоритма [90] и представляющий собой последовательность локальных модификаций сетки для достижения ее квази-равномерности.
Нужно отметить, что теоретического обоснования, почему сходится адаптивный алгоритм построения квази-оптимальных сеток для функций с особенностями, нет. Однако в Теореме 2.18 автору удалось построить оценки, показывающие возможность уменьшения нормы ошибки при интерполяции гладких функций [11, 53]. Вопрос теоретического обоснования в общем случае не решен, так же как и вопрос апостериорных оценок ошибки в рамках данной технологии, хотя для случая изотропного сгущения предложен конструктивный метод построения таких оценок [63], а для некоторых краевых задач и некоторых типов анизотропных сеток предложены оценки [174]. Альтернативная адаптивная технология [260, 59, 161], основанная на локальном измельчении/огрублении сетки на основе выравнивания апостериорно оцененной ошибки, естественным образом интегрирует оценку ошибки, однако использует иерахические локально измельчающиеся сетки, элементы которых наследуют форму начальной сетки.
Общий адаптивный алгоритм весьма технологичен, поскольку его структура предполагает выполнение четырех совершенно независящих друг от друга шагов:
шаг инициализации: построение начальной сетки;
нахождение сеточного решения (непрерывной кусочно-линейной функции);
восстановление сеточного гессиана;
генерация сеток, квази-равномерных в заданной метрике.
Ключевым свойством алгоритма является полная автономность от пользователя шагов 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] решения смешанных конечно-элементных систем:
на входе дается линейная смешанная конечно-элементная система в дисассемблированном (по ячейкам сетки) виде, а также координаты вершин симплексов и граф связности вершин;
с использованием геометрических данных система преобразуется локально по ячейкам в гибридную систему за счет добавления множителей Лагранжа на грани симплексов;
в расширенной системе локально по ячейкам исключаются скалярные неизвестные внутри симплексов и потоки на гранях симплексов, порождая симметричную положительно определенную разреженную матрицу;
итерационно решается разреженная конденсированная система для множителей Лагранжа;
восстанавливаются локально по ячейкам скалярные неизвестные внутри симплексов и потоки на гранях симплексов.
Ключевым наблюдением является то, что все этапы, кроме четвертого, арифметически
дешевы и легко параллелезуемы, в силу локальности действий. Параллельное решение
линейной системы, как наиболее трудоемкий этап всей цепочки, осуществляется быстрос-
ходящимся многосеточным методом [102, 103, 16], реализация которого позволила сфор
мировать всю технологическую цепочку.
^ Таким образом, алгоритмы, изложенные в третьей главе, открывают путь к исиоль-
f зованию параллельных адаптивных технологий для некоторых практически значимых
неконформных аппроксимаций.
Алгоритмы четвертой главы объединяет применение и тестирование в приложениях на прямоугольных сетках, хотя использование ни одного из методов не ограничено классом прямоугольных сеток. Рассматриваются две основные задачи гидродинамики: уравнения Навье-Стокса и задача многофазной фильтрации. Являясь нелинейными нестационарными уравнениями, они допускают ряд технологий получения решения [216, 132, 156, 243, 198, 145, 3, 26, 4]. Здесь мы ограничимся двумя наиболее востребованными методами. Во-первых, это полностью неявные аппроксимации по времени, обеспечивающие абсолютную устойчивость к шагу по времени и порождающие на каждом временном шаге систему нелинейных уравнений, которую необходимо решить. Во-вторых, это проекционный алгоритм для уравнений Навье-Стокса, сводящийся к решению одной или нескольких линейных систем, с нежесткими ограничениями на шаг по времени [198, 154, 248]. Преимуществом первого метода является возможность управлять шагом по времени на основе внешних данных (правой части, граничных условий), а не внутренней динамики системы, что весьма привлекательно в инженерных приложениях. Преимущества второго метода — в упрощении решаемых систем уравнений и существенном ускорении их решения на каждом шаге, что обеспечивает этим методам наибольшую скорость расчетов [248].
Проекционные алгоритмы были предложены в работах [101, 242, 42] и их смысл сводился к двухшаговой процедуре: с помощью уравнения моментов предсказать поле скоростей на очередном шаге по времени и спроектировать это поле в пространство соле-ноидальных функций. Предсказание возможно за счет явной аппроксимации по времени уравнения моментов [101] (при этом налагаются существенные ограничения на временной шаг), за счет метода характеристик [216] (параллельная реализация которого очень сложна [45]), за счет полунеявной аппроксимации по времени уравнения моментов [198, 170], чье единственное отличие от неявной аппроксимации — линеаризация нелинейного конвективного члена вида (ufc+1 V)ufe+1 —> (ик V)ufc+1. Последний подход предполагает решение линеаризованного уравнения моментов, представляющего собой в трехмерном случае три сингулярно возмущенных уравнения реакции-конвекции-диффузии для каждой из компонент скорости. Именно это приложение потребовало разработки эффективного параллельного метода решения сингулярно возмущенного уравнения конвекции-диффузии.
Сингулярно возмущенные операторы предоставляют дополнительные возможности для эффективной параллелизации решателей на основе методов декомпозиции с перекрывающимися подобластями. В применении к параболическим уравнениям наличие сингулярно возмущенного члена нулевого порядка привело к ряду родственных методов [176, 177, 74], основанных на экспоненциальном убывании сеточной функции Грина.
Предлагаемый алгоритм базируется на важном свойстве сингулярно возмущенного уравнения конвекции-диффузии: возмущения решения, вызванные локальным возмущением краевого условия, распространяются анизотропно. В случае постоянного вектора переноса точечное возмущение распространяется только вниз по потоку, быстро убывая в поперечном направлении и вверх по потоку. Это свойство может быть переформулировано в терминах анизотропного поведения сеточной функции Грина [165, 211, 137], которое может быть использовано при построении как декомпозиционных [222, 137], так и многосеточных методов [35, 212].
На основе этих свойств оператора конвекции-диффузии конструируется параллельный двухуровневый метод Шварца [137, 141], быстро сходящийся даже для малых налеганий между подобластями. На первом уровне область разбивается на налегающие полосы (слои), перпендикулярные доминирующему направлению переноса, и определяются итерации Шварца в направлении вниз по потоку. На втором уровне каждая подобласть-полоса разбивается на налегающие подобласти-квадраты и задаются несколько итераций Шварца между квадратами с четными и нечетными индексами. Естественный параллелизм методу дает именно второй уровень, поскольку каждая подобласть-квадрат может обрабатываться своим процессором. Если вектор переноса сильно отклоняется от заданного направления (например, вследствие завихрений), разбиение первого уровня должно быть адаптировано к поведению вектора.
Рассмотренный метод является первым итерационным методом для оператора конвекции-диффузии с плоско-параллельным полем переноса, чья скорость сходимости теоретически ограничена сверху константой, не зависящей от малости коэффициента диффузии е и шага сетки h, если є < h, см. Теорему 4.1, впервые доказанную в [136]. Метод обобщен автором на случай полей переноса с локальными завихрешюстями, при этом в Теореме 4.6 автором доказано, что итерационная скорость сходимости остается не зависящей от малости є и h при разумных предположениях на разбиение на подобласти. Арифметическая цена одной итерации этого метода (как и его инициализации) пропорциональна числу неизвестных в системе, поскольку представляет собой набор линейных систем с матрицами ограниченного порядка (в двумерном случае) или ленточными матрицами с малой шириной ленты (в трехмерном случае). Более того, метод легко параллелезуем. Отметим, что предложенный метод, будучи первым проанализированным декомпозиционным алгоритмом, универсальным по отношению к сеточному числу Пекле, получил обобщения на нелинейные задачи [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 decomposition 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], который преобразует исходную систему к специальному блочному виду, стремясь максимально отщепить зависимость переменной "давление" от переменных "концентрации". Такой подход имеет два потенциальных преимущества: во-первых, можно строить хороший переобуславливатель только для самых жестких блоков, переобуславливая другие блоки вычислительно дешевыми технологиями. Во-вторых, этот подход не зависит от структуры сетки, поскольку пере-обуславливатели для блоков могут опираться только на их алгебраическую структуру.
Второй обсуждаемый подход — многосеточный метод с ускоренным разгрублением SCMG [267, 184], который использует иерархичность трехмерной прямоугольной сетки. Несмотря на структурное подобие SCMG и стандартного V-цикла геометрического многосеточного метода MG [40, 158], они опираются на разные итерационные механизмы. Трехмерный метод MG использует последовательные коррекции на более грубых уровнях, в то время как SCMG — переобуславливатель без пост-сглаживания с коррекцией (аналогичной методу агрегирования из главы 1) в подпространстве двумерных функций, где используется усиленный (со многими сглаживаниями, в силу малости задач) двумерный метод MG для решения агрегированной подзадачи. Именно это обстоятельство позволяет использовать в SCMG простейшие геометрические операторы интерполяции и сужения даже для приложений с сильно меняющимися коэффициентами, что является неприемлемым для классического метода трехмерного метода MG [115]. Несмотря на разные механизмы, лежащие в основе используемых методов (возможность отщепления давления и эффективность вертикального агрегирования), численные эксперименты показали эффективность обоих методов решения возникающих линейных систем.
Перейдем к краткому содержанию диссертационной работы, состоящей из четырех глав.
В первой главе, состоящей из четырех разделов, рассматриваются некоторые блочные технологии решения неструктурированных систем. Некоторые известные методы представлены в рамках обзора; некоторые новые технологии представлены без численного анализа алгоритмов, но с экспериментальными данными; метод агрегирования рассмотрен подробно как в теоретическом, так и в экспериментальном плане.
Обзор современных последовательных методов решения неструктурированных систем представлен в разделе 1.1. Рассмотрены основные асимптотические характеристики методов точной факторизации, неполной факторизации, алгебраических многосеточных методов и метода фиктивного пространства. Характерные черты современных реализаций методов проиллюстрированы на последовательностях матриц возрастающей размерности и заданной структуры разреженности.
Далее рассматриваются новые комбинации известных технологий, оказывающиеся весьма эффективными для неструктурированных систем специального вида. Быстрое приближенное решение систем с матрицами тензорного ранга 2, порождаемых неструктурированными призматическими сетками, рассмотрено в разделе 1.2. Предлагаемый метод на практике характеризуется высокой скоростью сходимости и квази-оптимальной арифметической ценой.
Технология переупорядочивания неизвестных, позволяющая для задач с анизотропными коэффициентами эффективно использовать легко параллелезуемый блочно-диагональный переобуславливатель, излагается в разделе 1.3.
В разделе 1.4 рассмотрен и проанализирован метод агрегирования как один из декомпозиционных алгоритмов, позволяющий использовать в качестве блоков последовательные методы решения (любые из рассмотренных выше), обладающий очень удобной структурой для параллелизации и сконструированный для эффективного параллельного решения задач на неструктурированных сетках и/или с гетерогенными коэффициентами. Там же рассмотрены некоторые приложения этого метода.
Результаты первой главы сводятся к следующему: экспериментально исследована
асимптотика арифметической работы известных технологий решения неструктуриро
ванных систем; представлен новый метод эффективного решения систем с матрицами
тензорного ранга 2; рассмотрена и протестирована технология формирования блочных
переобуславливателей для неструктурированных задач с анизотропным коэффициентом
диффузии; предложен и исследован теоретически и экспериментально метод агрегиро-
* вания, удобный для параллелизации и применимый для задач с гетерогенными коэффи-
циентами, аппроксимированных на неструктурированных сетках.
Во второй главе, состоящей из пяти разделов, рассматриваются свойства оптимальных симплициальных (треугольных и тетраэдральных) сеток, а также свойства и методы построения их аппроксимаций, квази-оптимальных сеток.
В разделе 2.1 представлен теоретический анализ этих сеток: введены основные понятия, доказано существование оптимальных сеток в двумерном случае, для норм ошибки показаны оценки снизу на оптимальных сетках как в норме Loo > так и в норме Lp, р Є (0, +оо], и оценки сверху на квази-оптимальных сетках в норме Ь^.
В разделе 2.2 рассмотрены основные технологии адаптивного алгоритма построения квази-оптимальных сеток: представлены способы восстановления сеточного гессиана и их анализ, методика построения сеток, квази-равномерных в заданной метрике, а также приведен сам адаптивный алгоритм. Для случая минимизации ошибки кусочно-линейной интерполяции гладкой функции рассмотрен механизм сходимости адаптивного алгоритма.
Одно из основных преимуществ представленного алгоритма, простота управления адаптацией за счет модификации метрики, рассмотрено в разделе 2.3. Выведены оценки сравнения ошибок на сетках, квази-равномерных в модифицированной метрике, и ошибок на квази-оптимальных сетках.
В разделе 2.4 обсуждаются сложности в адаптации, порождаемые дискретным заданием границы расчетной области, и предлагается и анализируется метод квадратичного восполнения кусочно-линейной поверхности, частично снимающий остроту проблемы.
В разделе 2.5 излагается параллельный метод адаптивного построения квази-оптимальных сеток, включающий как параллельную генерацию сеток, так и параллельное решение сеточных задач, а также рассматривается влияние управления адаптацией на параллельное ускорение.
Результаты второй главы сводятся к следующему: представлен теоретический анализ оптимальных симплициальных сеток (существование и асимптотические свойства ошибки интерполяции); теоретически исследованы свойства их приближений (квазиоптимальных сеток); представлен и проанализирован новый способ восстановления гессиана сеточной функции; предложен теоретический базис для вывода оценок ошибки на сетках с управляемой адаптацией; предложен, проанализирован и протестирован новый метод квадратичного восполнения кусочно-линейной поверхности; рассмотрен и проанализирован параллельный метод адаптивного построения квази-оптимальных сеток, и исследовано влияние управления метрикой на параллельное ускорение.
В третьей главе, состоящей из пяти разделов, рассматриваются блочные параллельные технологии решения систем линейных уравнений, возникающих в некоторых прикладных неконформных аппроксимациях.
В разделе 3.1 дается макро-гибридная формулировка краевых (диффузионных) задач на нестыкующихся сетках, а также общая технология параллельного итерационного решения соответствующих линейных систем.
В разделе 3.2 рассматриваются интерфейсные переобуславливатели, являющиеся ключевыми аспектами этой технологии. Строятся два типа переобуславливателей, применимых к широкому классу конформных сеток в подобластях, и универсальных по отношению к скачку коэффициента диффузии, количеству подобластей и малости коэффициента реакции.
Сочетание параллельной итерационной технологии с адаптивной генерацией расчетных сеток по подобластям независимо друг от друга порождает новый класс параллельных адаптивных алгоритмов, описанный в разделе 3.3.
В последних разделах главы рассматриваются неконформные аппроксимации на конформных сетках. В разделе 3.4 предлагается и исследуется мозаично-скелетонный пере-обуславливатель в декомпозиционном алгоритме для макро-гибридной формулировки с неконформными множителями Лагранжа.
В разделе 3.5 рассмотрен параллельный многосеточный алгоритм для решения неконформных конечно-элементных систем, эквивалентных компактной алгебраически конденсированной форме записи смешанных конечно-элементных систем для эллиптических уравнений второго порядка.
Результаты третьей главы сводятся к следующему: технология с внутренним итерационным процессом для параллельного итерационного решения систем, порождаемых макро-гибридными формулировками на нестыкующихся сетках, обобщена и обоснована на случай произвольных конформных симплициальных сеток; предложена и исследована новая декомпозиционная технология параллельного итерационного решения таких систем; предложена и протестирована новая параллельная технология адаптивного решения краевых задач на неструктурированных сетках; рассмотрен новый метод декомпозиции с использованием мозаично-скелетонного аппроксимации плотных матриц в качестве интерфейсного переобуславливателя; предложена и протестирована новая параллельная технология решения смешанных конечно-элементных систем, порождаемых иерархическими локально сгущающимися сетками.
В четвертой главе, состоящей из четырех разделов, рассматриваются блочные параллельные технологии решения систем линейных уравнений, порождаемых аппроксимациями на трехмерных прямоугольных сетках.
В разделе 4.1 формулируется и анализируется параллельный двухуровневый метод Шварца в применении к сингулярно возмущенному уравнению конвекции-диффузии.
Различные параллельные приложения этого метода в нестационарных задачах вы-
числительной гидродинамики рассматриваются в разделе 4.2. Наряду с уравнениями конвекции-диффузии, эти задачи порождают еще один класс проблем: необходимость решать последовательность больших жестких систем, возникающих на каждом шаге по времени.
В разделе 4.3 представлено несколько адаптивных итерационных технологий, ускоряющих время решения последовательности систем, порождаемых в процессе аппроксимации по времени. Эти технологии легко параллелезуются и не используют прямоуголь-ность расчетной сетки; их уместность в этой главе диктуется исключительно соображениями удобства читателя, только что познакомившегося с основными постановками задач и алгоритмами их сведения к набору модельных подзадач.
В разделе 4.4 излагаются параллельные технологии эффективного решения систем уравнений многофазной фильтрации, как многоуровневого, так и блочного типов. Во всех случаях фактором, определяющим высокую скорость итерационной сходимости, является использование свойств операторов задачи.
Результаты четвертой главы сводятся к следующему: параллельный двухуровневый метод Шварца в применении к сингулярно возмущенному уравнению конвекции-диффузии обобщается и обосновывается на случай произвольных полей вектора скорости с доминирующим направлением конвекции; рассмотрены численные параллельные трехмерные приложения разработанной технологии (нестационарные уравнения Навье-Стокса, уравнения тепло-массопереноса в простейшем химическом реакторе); рассмотрены и протестированы три параллельные стратегии адаптивного выбора начального приближения при итерационном решении последовательностей линейных систем с одной матрицей, разными матрицами, и нелинейных систем; разработаны и исследованы две параллельные технологии решения линейных систем, возникающих при полностью неявных аппроксимациях уравнений многофазной фильтрации.
Основным научным результатом диссертации является разработка новых эффективных технологий параллельного решения краевых задач с использованием адаптивных сеток. Автор выносит на защиту следующие научные результаты:
предложен, проанализирован и реализован новый параллельный метод агрегирования, удобный в параллельной реализации, доказана его универсальность по отношению к гетерогенности коэффициента диффузии и количеству подобластей;
представлен теоретический анализ оптимальных симплициальных сеток; теоретически исследованы свойства их приближений, квази-оптимальных сеток; доказаны
* оценки сверху и снизу для интерполяционных ошибок на этих сетках;
3. реализованы и расширены основные компоненты технологии параллельного адап-
тивного алгоритма построения квази-оптимальных сеток на основе восстановления гессиана (новый метод восстановления гессиана сеточной функции с особенностями и его теоретический анализ, последовательные и параллельные методики построения сеток, квази-равномерных в заданной метрике, тестирование технологии на ряде приложений);
рассмотрены теоретические и практические вопросы управления адаптацией в рамках вышеупомянутой адаптивной технологии; доказаны оценки сверху для интерполяционных ошибок на сетках с управляемой адаптацией;
предложены, обоснованы и реализованы два параллельных алгоритма решения систем, порожденных макро-гибридными формулировками на нестыкующихся сетках; доказана их универсальность по отношению к количеству подобластей; на основе этих алгоритмов предложена и реализована новая параллельная технология параллельного адаптивного решения краевых задач;
предложен, обоснован и реализован параллельный двухуровневый метод Шварца для решения сингулярно возмущенного уравнения конвекции-диффузии с доминирующим направлением поля скорости; доказана универсальность метода по отношению к большим числам Пекле; рассмотрен ряд его трехмерных приложений метода для решения нестационарных уравнений Навье-Стокса, уравнения тешю-массопереноса в простейшем химическом реакторе;
предложены, реализованы и протестированы параллельные технологии адаптивного выбора начального приближения при итерационном решении последовательностей линейных систем с одной матрицей, разными матрицами, и нелинейных систем;
разработаны и исследованы параллельные технологии решения линейных систем, возникающих при полностью неявных аппроксимациях уравнений многофазной фильтрации.
Практическая ценность разработанных адаптивных технологий состоит в их применимости для решения широкого круга краевых задач математической физики. Важным свойством технологий также является модульность технологических составляющих, что позволяет использовать их как в уже реализованных технологических цепочках в качестве ингредиента, так и создавать новые технологические цепочки. Актуальность парал-лелизации технологий заключается в их приложении к трехмерным задачам, требующим большого числа степеней свободы.
Основные результаты диссертации докладывались автором: на всероссийской школе-конференции (Казань, 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 Petrol, 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, И, 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 базируется на экспериментах автора с неопубликованным методом, предложенным и реализованным С.А.Горейновым, исходя из постановки проблемы автором.
Автор искренне признателен всем соавторам и коллегам за сотрудничество. Автор сердечно благодарит Ю.А.Кузнецова за многолетнюю поддержку, и С.В.Непомнящих за
^
плодотворное обсуждение многих научных вопросов.
Блочные переобуславливатели и переупорядочивание для задач с анизотропными коэффициентами
Ключевым свойством алгоритма является полная автономность от пользователя шагов 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 обоснована независимость скорости сходимости итерационных алгоритмов от шага сетки (методы второй группы), а также числа подобластей и величины коэффициентов диффузии и реакции для уравнений реакции-диффузии (методы обеих групп). Построены параллельные адаптивные алгоритмы и на численных примерах проиллюстрированы их основные свойства.
Параллельная генерация сеток, квази-равномерных в заданной метрике
Для графов, чьи степени вершин ограничены снизу и сверху, выбор а управляет минимально возможной степенью внутри блока. Выбор а — 0 снимает это ограничение.
Выбор /3 задает тип шаблона в блоках в следующем смысле. Если шаблон сеточного оператора является семиточечным (конечные разности, или конденсированные матрицы для смешанных конечных элементов), то выбор /3 2/3 порождает блоки с семиточеным внутренним шаблоном; выбор 1/3 /3 2/3 допускает внутренние шаблоны с не менее, чем 5 точками; выбор /3 1/3 позволяет иметь внутренние шаблоны с не менее, чем 3 точками. Выбор /3 = 0 снимает это ограничение.
Выбор 7 управляет влиянием матрицы весов на создание блоков, т.е. действительная разреженность блока зависит от матрицы весов. Например, если элементы матрицы весов W имеют только два значения, w и w , w » го 0, и ненулевая структура блока должна быть определена г//-ребрами, то достаточно выбрать 7 w /w .
Предположим, что исходная матрица А после переупорядочивания и разбиения на блоки Жадным алгоритмом имеет блоки {Bkj}bkj=1. Будем использовать переобуславли-ватель на основе блочного метода Гаусса-Зейделя Разумеется, возможны и другие типы переобуславливателей. Например, блочный метод Якоби, blockdiag{5jj}, является легко параллелезуемым, хотя и имеет более низкую скорость сходимости. Альтернативой могут быть подходы на основе неполной блочной факторизации, впервые рассмотренной в [202]. Кроме того, точная факторизация блоков ВЦ может быть заменена неполной факторизацией.
Областью применения предложенного алгоритма являются матрицы, порожденные сеточными операторами с анизотропией. Источником анизотропии может быть либо анизотропия коэффициентов диффузии в операторе диффузии, либо анизотропия ячеек сетки, либо другие причины. В этом разделе мы рассмотрим случай анизотропных коэффициентов диффузии для смешанной конечно-элементной аппроксимации оператора с краевыми условиями Дирихле в прямоугольных двумерных и трехмерных областях на неструктурированных локально сгущающихся симплициальных сетках (Рис.1.4, 1.5). Тензор диффузии диагоналей, и коэффициенты могут терпеть скачки на границе между двумя равными подобластями: К = dia,g{kxx, куу, kzl} в трехмерном случае, K=diag{fcl куу} в двумерном случае, і — 1,2. Рассматриваемые матрицы порождаются алгебраической конденсацией смешанной конечно-элементной матрицы. Для неструктурированных сеток структура разреженности совершенно нерегулярна, хотя в каждой строке матрицы не более 7 ненулевых элементов, а в двумерном случае - не более 5 ненулевых элементов.
Во всех экспериментах весовая матрица содержит модули внедиагональных элементов конденсированной матрицы. Параметры ограничений выбраны как а = 0.9, (3 — 0.1, 7 — 0.4 (2D), 7 = 0.1 (3D). В таблицах представлены характеристики итерационнно-го решения GMRES(20) конденсированных систем со случайными правыми частями, с подавлением начальной невязки в 104 раз, с переобуславливателем и без: время инициализации переобуславливателя ипц, и время итераций Utr. Временные характеристики получены на персональном компьютере PC Intel Pentium 4, 2.5 ГГц.
Рассмотрим прямоугольник [0,2] х [0,1] и регулярную триангуляцию, сгущающуюся к границе между подобластями х — 1 (Рис. 1.4). Отношение максимального и минимального размера треугольников около 1000. Алгебраически конденсированная матрица СМКЭ аппроксимации оператора (1.9) имеет порядок 12330 и не более 5 ненулевых элементов в строке. Время работы Жадного алгоритма не превышает 0.01 секунды во всех случаях его применения. Результаты Таблицы 1.5 показывают, что число итераций непереобу-словленного GMRES очень велико, что обусловлено большой жесткостью матрицы из-за малого шага сгущающейся сетки. Кроме того, в зависимости от свойств коэффициентов число блоков #block, порождаемых Жадным алгоритмом, сильно разнится. Большие блоки появляются во всех случаях, и их максимальные размеры 5116 (случаи 6,7), 9208 (случаи 3,4) и 12258 (случаи 1,2,5). Однако, только случаи 1,2,5 имеют чрезмерно большое время инициализации, хотя и обеспечивают очень хорошую скорость сходимости. Причина такого расхождения в том, что в случаях 1,2,5 блоки имеют около 5 ненулевых элементов в строке, а в остальных случаях - не более трех. Таким образом, на эффективность инициализации влияет не размер блока, а степень его разреженности. Наконец, уместно отметить положительный эффект применения блочно-диагонального переобу-славливателя в случае сильно анизотропных коэффициентов диффузии.
Рассмотрим параллелепипед [0,2] х [0,1] х [0,1] и регулярные тетраэдризации, сгущающиеся к границе между подобластями х = 1 (Рис.1.5). Отношение максимального и минимального размера тетраэдров около 250. Алгебраически конденсированная матрицаСМКЭ аппроксимации оператора (1.9) имеет порядок 40254 и не более 7 ненулевых элементов в строке. Время работы Жадного алгоритма не превышает 0.04 секунды во всех случаях его применения.
Результаты численных экспериментов приведены в Таблице 1.6. Аналогично двумерному случаю, число итераций непереобусловленного GMRES тоже очень велико, что также объясняется малостью шага сгущающейся сетки. Однако, в отличие от двумерного примера, количество блоков, порожденных Жадным алгоритмом, не сильно зависит от анизотропии коэффициентов: средний размер блоков варьируется в пределах от 5 до 10. Более того, большие блоки не производятся вовсе: для почти изотропных случаев, самые большие блоки имеют порядок 300-400 с 4 ненулевыми элементами в строке, а в остальных случаях их размер не превышает 252 и они имеют такую же разреженность. Поэтому итерационное решение становится весьма эффективным даже в случаях слабой анизотропии. Отметим существенное ускорение итерационного решения в терминах как числа итераций (более 100 раз), так и времени итерационного решения (более 50 раз).
Отметим, что хотя Жадный алгоритм является последовательным, он очень экономичен как с точки зрения арифметической работы, так и с точки зрения запросов на компьютерную память. Поэтому его можно применять для обработки очень больших матриц на одном процессоре. В свою очередь процесс итерационного решения может быть распараллелен путем замены блочного метода Гаусса-Зейделя на блочный метод Якоби. При этом блоки могут быть распределены по процессорам с учетом арифметической работы для инициализации и применения переобуславливателя. Будучи менее "связывающим".
Переобуславливатели с внутренним итерационным процессом для множителей Лагранжа
Мы рекомендуем рассматриваемую в [234] технику расцвечивания и переупорядочивания по следующей причине. Матрица А разрежена за исключением т строк, содержащих много элементов. Следовательно, она может быть эффективно разбита на блоки жадным алгоритмом расцвечивания. На практике число получаемых блоков не превышает десяти. Расцвечивание порождает такое переупорядочивайте строк и столбцов А, что диагональный блок, соответствующий какому-то цвету — диагональная матрица, и как последовательная, так и параллельная реализации BSOR очень просты. Число итераций BSOR является нечувствительным к скачкам коэффициента р(х). Отметим, что все ингредиенты метода очень удобны для параллельной реализации и приложений на неструктурированных сетках.
Прежде всего, сформулируем предположения, в рамках которых параллелизация итерационного решения методом агрегирования может быть эффективной.
Во-первых, пусть мы не можем держать глобальную матрицу жесткости и строить переобуславливатель для нее на одном процессоре (сетка уже с миллионом тетраэдров требует около 100 мегабайт для матрицы жесткости и около 400 мегабайт для структуры данных алгебраического многосеточного метода), однако предполагается, что матрица на грубой сетке А может храниться на одном процессоре (для нашей сетки и разбиения на две подобласти она займет не более 4 мегабайт).
Во-вторых, мы предположим, что число процессоров в нашем компьютере не очень велико, поскольку иначе порядок агрегированной матрицы будет приближаться к порядку глобальной матрицы, что запрещено предположением 1.
В-третьих, предположим, что всякий итерируемый вектор v Є R" разбит на т примерно равных непересекающихся подмножеств; г-е подмножество и г-я блочная строка [Ал,... ,Агт] обрабатываются только одним процессором с индексом г. Умножение матрицы на вектор требует векторов Vj, которые не известны на процессоре с индексом г. Эти данные обеспечиваются общепринятой методикой фиктивных узлов, которая присоединяет к г-му процессору копии тех элементов векторов Vj, і ф ];, которые необходимы для вычисления AijVj. Перед каждым умножением на матрицу, значения в фиктивных узлах должны быть обновлены соответствующими значениями на смежных процессорах (процедура Update). Другими словами, Update — единственная процедура, используемая при межпроцессорных обменах. Она может быть эффективно реализована через, например, асинхронные send/receive МРІ-функции.
В-четвертых, предположим, что арифметическая цена применения переобуславлива-телей в подобластях Вц предполагается монотонной функцией от щ. Это означает необходимость равномерного распределения степеней свободы по имеющимся процессорам для достижения их равномерной загрузки.
Реализация общепринятых вариационных методов итерационного решения линейных систем требует трех параллельных операций: умножение на матрицу, применение пере-обуславливателя, и вычисление скалярного произведения. Применение иереобуславли-вателя подразумевает три матрично-векторных произведения с матрицами В±, А и В. Матрица В\ блочно-диагональна, поэтому прозведение ее на вектор не требует никаких межпроцессорных обменов. Матрица В — результат нескольких BSOR итераций для агрегированной системы. Реализация BSOR итераций предполагает, что агрегированная матрица создана и разбита (после переупорядочивания) на блоки, причем диагональные блоки являются диагональными матрицами. Эти операции выполняются на этапе инициализации: процессор с индексом і вычисляет локальную матрицу агрегирования % и (щ +1) х п блочную строку матрицы А, соответствующую г-й подсетке. Блочная строка может быть ассемблирована на процессоре с индексом root для формирования матрицы А. Матрица А — разреженная, за исключением га строк, в которых много ненулевых элементов. Поэтому она может быть эффективно разбита на блоки стандартным последовательным жадным алгоритмом расцвечивания [234]. Как уже было отмечено, эта техника фактически производит такое переупорядочивание матрицы А, что любой диагональный блок, соответствующий какому-либо цвету — диагональная матрица. Рассылка цветов неявно определяет разбиение на блоки матрицы А для каждого процессора. Известно, что с оптимальным параметром релаксации и (вычисленным адаптивно на этапе инициализации [168, 272]), BSOR итерации являются эффективным сглаживате-лем [64, 272]. Кроме того, параллельные итерации BSOR используют только процедуру Update для межпроцессорных коммуникаций [64].
Отметим, что использование коррекции блочно-диагонального переобуславливателя за счет приближенного решения агрегированной системы играет троякую роль. Во-первых, мы можем применить очень простую последовательную процедуру расцвечивания для агрегированной матрицы, поскольку ее размер и число элементов существенно меньше, чем в исходной глобальной матрице. Во-вторых, по этой же причине существенно снижаются затраты на сглаживание BSOR итерациями. В-третьих, для жестких систем, жесткость агрегированной матрицы существенно меньше, чем у исходной матрицы, поэтому BSOR итерации оказываются эффективным сглаживателем. В качестве начального этапа численного тестирования метода (1.17) сравним его характеристики со стандартным аддитивным методом Шварца. Рассмотрим квази-равномерную неструктурированную тетраэдризацию единичного куба, состоящую из Л д = 1744 тетраэдров. Для измельчения сетки будем использовать равномерное разбиение каждого тетраэдра на 8 подтетраэдров. В задаче (1.10) положим р(х) — 1 и f(x) — 1. Разбиение на две подобласти с минимальным наложением 6 — h обеспечивается разделением узлов сетки плоскостью х\ — 0.5. В Таблице 1.7 представлены основные характеристики гибридного метода Шварца В и аддитивного переобуславливателя метода Шварца В\ в предположении, что В = Л-1, ВЦ = АЦ, І — 1,2. Фактически, матрица грубого уровня А не обращается, точное решение системы уравнений моделируется большим числом (30) BSOR итераций. Сходимость измеряется числом итераций метода сопряженных градиентов, необходимых для уменьшения начальной невязки в 106 раз. Число итераций с переобуславливателем В в 2.5 раза меньше, чем с переобуславливателем В\. Однако линейная зависимость числа обусловленности от 5 1 отмечается в обоих случаях.
Во втором эксперименте мы покажем, что наличие сглаживателя В в (1.17) устраняет (или уменьшает) чувствительность метода к числу подобластей. Рассмотрим тетраэдризацию кубической сетки в единичном кубе U — (0,1)3, h = 2 6, NT — 1572864. Опять обратимся к идеализированной ситуации: Вц — Ац, В — A l. Наряду с задачей (1.10).
Проекционный алгоритм для уравнения тепло-массопереноса в простейшем химическом реакторе
В качестве иллюстрации рассмотрим тестовый пример из [13, 189, 187] и объясним результаты численных экспериментов с теоретической точки зрения. Пусть и — решение первой краевой задачи для уравнения Пуассона с однородным краевым условием в области Г2 со входящим трехгранным углом, 17 = (О, І)3 \ [0,0.5]3: где f(x) —правая часть с особенностью, f(x) — 1/\х — х0\, и х0 = (0.5,0.5,0.5). Свойства решения и исследованы в [57]. Решение обладает слабой анизотропной особенностью около ребер трехгранного угла и сильной изотропной особенностью во входящем угле Хо из-за особенности правой части.
Рассмотрим три различные метрики. Первая метрика порождается гессианом, восстановленным из дискретного решения. Вторая и третья метрики получаются из первой с помощью методов управления (2.64) и (2.65), соответственно. Пусть весовая функция равна и(х) — l/f(x). Весовая функция усиливает адаптацию сетки вдоль ребер, поскольку влияние сильной точечной особенности в окрестности Хо не учитывается в весовой метрике. В этом случае мы измеряем ошибку в области D — П \ В, которая не содержит шар В радиуса го = 0.1 с центром в XQ.
Как показано на Рис. 2.5, ошибки для всех трех метрик пропорциональны J\f(Qh) 2 , что является асимптотически оптимальным результатом. Ошибка на квази-равномерных сетках, порожденных равномерной (единичной) метрикой, не оптимальны. След квазиоптимальной сетки показан на Рис. 2.6 слева. Для второй (изотропной) метрики, наихудший тетраэдр Д в ІЯ І-квази-равномерной сетке расположен вблизи одного из ребер трех входящих двугранных углов, см. Рис. 2.6 в центре. Однако, поскольку анизотропная особенность вдоль ребер слаба, множитель а(Н,Н) в Теореме 2.19 достаточно мал, и асимптотическое убывание ошибки аналогично случаю квази-оптимальных сеток. В свою очередь, взвешенный гессиан (2.65) сохраняет основные свойства исходного гессиана везде, кроме окрестности точечной особенности во входящем угле, куда сетка более не сгущается адаптивно, см. Рис. 2.6 справа. В этом случае для области D множитель а(Н, Н) близок к единице, и ошибка в области D ближе к ошибке на сетках, построенных на основе немодифицированного гессиана.
Во многих прикладных расчетах точная параметризация криволинейных поверхностей, представляющих части границы, неизвестны. Эти поверхности описываются треугольными сетками в пространстве (например, сетки, порождаемые системами САПР), что ухудшает работу адаптивных методов из-за недостаточного разрешения поверхности границы. Одно из возможных решений — использовать результат адаптивного расчета в качестве обратной связи для модели САПР. Подобный подход требует прямого вмешательства пользователя и может быть слишком сложным для некоторых приложений. Альтернативный подход основан на том обстоятельстве, что если поверхности границы достаточно гладкие (или кусочно-гладкие), то исходная треугольная сетка несет дополнительную информацию об этих поверхностях. В этом разделе мы рассмотрим влияние неточности в представлении границы на работу адаптивной технологии и предложим новый метод восстановления поверхности, уменьшающий это влияние как с численной, так и с теоретической точек зрения.
Напомним, что основой нашей адаптивной технологии является построение симплици-альной сетки, квази-равномерной в некоторой заданной метрике. Генерация сетки представляет собой последовательность локальных модификаций элементов сетки, повышающих их качество, являющееся мерой квази-равномерности. На модификации должны налагаться определенные ограничения. Так, например, симплициальная сетка не может вывертываться (элементы не могут налегать друг на друга). Это достигается контролем над сохранением знака ориентированных объемов окружающих симплексов. Кроме того, если узел лежит (или появляется) на граничной поверхности, его движение должно быть ограничено этой поверхностью. В случае параметрического задания поверхности, это требование обеспечивается за счет привязывания степеней свободы узлов к параметрическим переменным. В случае дискретного (кусочно-линейного) задания поверхности, движение узла по поверхности реализуется по специальной технологии [123], описание которой не входит в рамки диссертации. Таким образом, наличие топологических ограничений не является препятствием для реализации адаптивного процесса, даже в отсутствие явной и точной параметризации криволинейной поверхности.
Однако неточность в представлении криволинейной поверхности существенно влияет па уменьшение ошибки в ходе адаптивных итераций, как будет показано ниже на численном примере. Точность аппроксимации кусочно-гладкой поверхности можно существенно повысить за счет восстановления поверхности более высокого порядка, чем заданная кусочно-линейная дискретизация. Существует несколько методов восстановления поверхности более высокого порядка на основе кусочно-линейных поверхностей (см. [142, 200, 203, 204] и ссылки в них). В работе [200, 204] параметризация поверхности (а также другие характеристики поверхности) вычисляется из производных функций, задающих параметризацию. В работе [142, 200] дискретная кусочно-линейная поверхность аппроксимируется кусочно-квадратичной поверхностью методом наименьших квадратов. Метод, предлагаемый ниже, использует технологию дискретной дифференциальной геометрии для вычисления приближенного гессиана кусочно-квадратичной функции, задающей восстанавливаемую поверхность. Гессиан вычисляется на основе обобщенной формулировки, по аналогии с методом конечных элементов. Предлагаемый метод точен для квадратичных поверхностей.
Рассмотрим гладкий кусок поверхности Г, чья граница будет обозначаться через 9, и его кусочно-линейную аппроксимацию 1\ со своей границей вл.. Предположим, что узлы Г/і и Qh лежат на Г и Э, соответственно, хотя данное предположение не является необходимым на практике. Кусочно-квадратичная экстраполяция Г\ поверхностной триангуляции Г/( с гранями Tt определяется как непрерывная поверхность, являющаяся замыканием объединения открытых непересекающихся кусков Г\ локальных квадратичных экстраполяции над гранями Г(.
Локальная экстраполяция Г( описывается квадратичной функцией p2,t. Ниже мы будем опускать индекс t, если это не приводит к путанице. Опишем функцию ірч в удобной для нас локальной системе координат (i, 2), связанной с плоскостью грани Г(. В этой координатной системе многоточеная формула Тейлора для квадратичной функции ірі с гессианом Я 2 = {#$,2}piS=1 записывается как где fli, 02, Оз — вершины треугольника Tt, а РІ() — кусочно-линейная функция, такая что pi(aj) = Sij.
Восстановление гессиана Е4 2, основано, во-первых, на величинах a = (HV2ii,i), і = 1,2,3, представляющих проекции гессиана на ребра ii of Tt. Здесь и далее мы будем использовать как для сеточного ребра, так и для соответствующего вектора. В локальной системе координат, векторы задаются двумя координатами = (1\, 1\). Предположим, что вектор ІІ начинается в узле щ и заканчивается в щ+і, где а\ — а\. Тогда, по определению а имеем