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

Падий, Александр Викторович

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

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

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

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

1998 ГОД ЗАЩИТЫ

01.01.07 КОД ВАК РФ

Диссертация по математике на тему «Алгебраические многосеточные методы для задач с малым параметром»

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

российская академия наук р Г б {$3ститут вычислительной математики

2 1 ДЕК

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

ПАДИЙ Александр Викторович

АЛГЕБРАИЧЕСКИЕ МНОГО СЕТОЧНЫЕ МЕТОДЫ ДЛЯ ЗАДАЧ С МАЛЫМ ПАРАМЕТРОМ

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

автореферат

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

Москва - 1998

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

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

доктор физико-математических наук КУЗНЕЦОВ Ю. А. Официальные оппоненты:

доктор физико-математических наук ТЫРТЫШНИКОВ Е. Е. кандидат физико-математических наук КАПОРИН И. Е.

Ведущая организация:

Механико-математический факультет МГУ

Защита состоится "_"_199_г. в_часов на

заседании специализированного совета К.003.47.01 в Институте вычислительной математики РАН по адресу: Москва, ул. Губкина, 8.

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

Автореферат разослан "_"_199_г.

Ученый секретарь

специализированного совета

кандидат физико-математических наук С. А. ФИНОГЕНОВ

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы.

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

Цели работы.

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

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

3. Реализация разработанных методов на многопроцессорных ЭВМ.

Общая методика исследований.

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

Научная новизна.

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

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

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

Апробация работы.

Основные результаты диссертации докладывались на следующих конференциях и семинарах:

• на семинарах лаборатории вычислительной математики ИВМ РАН, семинарах Франко-Русского центра им .Ляпунова, семинарах в университете г. Ювяскюла (Финляндия) и в университете г. Наймегена (Голландия);

• на международной конференции "Algebraic Multilevel Iteration Methods", г.Наймеген, Голландия, 14-16 июня 1996 г., доклад " О стабилизированном аддитивном алгебраическом методе многоуровневых итераций";

• на международной конференции "Preconditioned Iterative Solution Methods for Large-Scale Problems in Scientific Computations", г.Наймеген, Голландия, 27-29 мая 1997 г., доклад "Об итерационном методе решения анизотропных эллиптических задач";

• на конференции Молодых ученых, Новосибирск, 20-22 апреля 1998 г., доклад "Модельный анализ приближенного варианта метода окаймления";

• на международном семинаре "Iterative Processes for Solving Equations", г.Киль, Германия, 3-5 июля 1998 г., доклад "Об эффективном параллельном методе решения задач теории упругости";

• на международной конференции "Mathematical Modelling and Computational Methods in Mechanics and Geodynamics", г.Прага, Чехия, 7-11 июля 1998 г., доклад "О многосеточном алгоритме решения задач теории упругости, скорость сходимости которого не зависит от регулярности используемых дискретизационных сеток";

• на международной конференции "Iterative Solution Methods for the Elasticity Equations as Arising in Mechanics and Biomecanics", г.Наймеген, Голландия, 28-30 сентября 1998 г., доклад "О параллельном многоуровневом методе для итерационного решения задач теории упругости";

Публикации.

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

Структура диссертации.

Диссертация состоит из введения, пяти глав и списка литературы (95 наименований). Общий объем работы: 80 страниц.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

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

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

Определим билинейную форму

A(u,v) = /ZajM^-^-, u,veV с Я'(Я),

u ij UXiUXj

где ii — многоугольник, а матрица коэффициентов [а^-(г)] является симметричной равномерно положительно определенной для всех х £ Q, и рассмотрим следующую вариационную задачу:

Найти такое и £ V С Я0' (Q), что

A(u,v) = (f,v), feLi(Q), Vu£V,

где под V = Vi понимается конечно-элементное подпространство (Q), принадлежащее последовательности пространств Vi С V-j С ... С С V/, заданных на последовательности триангуляций f2i С Пг С ... С iii-i С О/. Рассмотрим систему линейных алгебраических уравнений, порождаемую конечно-элементной дискретизацией исследуемой вариационной задачи:

A«x = f.

Для построения переобуславливателя к ASO воспользуемся подходом, основанным на рекурсивном устойчивом двухуровневом разбиении пространств сеточных функций. Под устойчивым двухуровневым разбиением пространства V будет пониматься такое разбиение V = V\ + V2, для которого существует такая постоянная 7 £ [0,1), что выполняется следующее усиленное неравенство Коши-Буняковского:

А{и, v) < 7{А(и, и) ■ A(v,v)y^ для всех и £ Vj, v £ V2.

Примером устойчивого друхуровневого разбиения V = + ^ является разбиение, получаемое в результате использования иерархического базиса конформных кусочно-линейных конечных элементов. Важной особенностью данного разбиения является то, что константа 7 может быть определена локально, на уровне суперэлементов сетки (из чего следует, что она не зависит от скачков в матрице коэффициентов при условии, что постоянна на каждом из элементов труангуляции самого грубого уровня измельчения) и отделена от единицы независимо от формы используемых треугольников. В случае использования для дискретизации прямоугольных сеток (в том числе неравномерных) 7 = 1/^2; а в случае использования триангуляции общего вида у = у/З/2.

Рассмотрим задачу, соответствующую двухуровневому разбиению:

} V,

} У2 '

Введем две вспомогательные матрицы Мц и М22, являющиеся аппроксимациями (переобуславливателями) к блокам Ац и А22 соответственно, и рассмотрим следующие способы построения переобуславливателя М к Л, основанные на блочном 2x2 разбиении матрицы:

• Мультипликативная схема:

МЫ1 =

• Аддитивная схема:

' Лп Л12 *

. Л21 л22 "2

'Мп 0 ' 1 Мп1Аи

Л21 М22 0 I

Мам =

Мп О О Мп

Теорема 1.1. Предположим, что выполняются следующие условия: у[ЛцУ1 < МцУ1 < (1+ А\\У\ для всех VI,

V2 Л2гУ2 < \\М^У2 < (1 + ¿г^ЛггУг для всех

где ¿1 и ¿2 — некоторые положительные постоянные. Тогда справедливы следующие соотношения:

с™;гVхЛV < уТМ,ппУ < СтцУтАу для всех V, СаА*УТЛУ < V7МыЫУ < С,ииУТАУ для всех V,

с константами

СппП — 1) Стц —

1+1 [¿1 + ¿2 + У(«5! - 52у +

Оий — , ) СаМ — . /, ч 1

1 + 7 Ло(1~7)

1 -72 До =

Рассмотрим построение аппроксимаций Мц и М22 к блокам Ац и А-щ. Как может быть легко показано, блок Ац хорошо обусловлен по отношению к шагу сетки: к(Ац) = 0(1), поэтому во многих случаях в качестве Мц может использоваться диагональ матрицы Ац: Мц = diag (Ац). Однако, данный подход неприменим, если используемая триангуляция содержит элементы, близкие к вырожденным, так как в этом случае обусловленность Ац резко ухудшается. Построение устойчивых к параметрам задачи переобуславлива-телей будет рассмотрено в Главах 2 и 4 настоящей диссертации. Поскольку блок А-22 соответствует матрице жесткости, порожденной оператором задачи при дискретизации на грубой сетке, в качестве М22 может быть выбран любой подходящий переобуславливатель. При рекурсивном использовании в качестве М-п двухуровневого переобуславливателя, аналогичного рассматриваемому, мы приходим к многоуровневой процедуре переобуславливания.

• Мультипликативная многоуровневая схема:

М®

м,ыь

I м'Г'ай1

Аддитивная многоуровневая схема:

М«

М®

А«

} V*

}

} •••

}

} V,

}

} V*-!

К сожалению, переобуславливание данного типа не приводит к итерационным методам, имеющим оптимальний порядок арифметической сложности, так как число обусловленности к (МСГ'а«) зависит от количества уровней измельчения. Для стабилизации числа обусловленности многоуровневого метода воспользуемся полиномиальной стабилизацией:

Мультипликативная схема с полиномиальной стабилизацией:

,(кГ1л(к)

м(1 =

М«

АЮ Л[к-1)

Л21 тпИ

I М£> АЪ О /

} Ук

где

Аддитивная схема с полиномиальной стабилизацией:

д/Ш _ 4(1) дЛ*) _

маЛ4 ~ л > таМ —

Л/}? о О

} и

} Ук-х

где

В качестве может выбираться либо чебышевский полином, наименее отклоняющийся от нуля на отрезке [а, Ь], где а и 6 — оценки границ спектра матрицы А/М либо многочлен с коэффициентами, порождаемыми ме-

тодом сопряженных градиентов.

Лемма 1.1. При использовании спектрально эквивалентных аппроксимаций Мух к блокам .а'п и пРи подходящем выборе степени полиномов щ (а также их коэффициентов) можно построить многоуровневый переобус-лавливателъ М« такой, что к ограничено независимо от I.

При этом оценка А^) не зависит от разрывов в матрице коэф-

фициентов [ау(х)]. Более того, в случае использования для дискретизации конформных кусочно-линейных конечных элементов оценка

«(Л/О-'лО) не

зависит и от параметров дискретизационной сетки.

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

} "новые" базисные функции } "старые" базисные функции '

А® 4? А®

Предлагается метод построения аппроксимаций Л/Ц , основанный на процедуре сохранения единственной (максимальной по модулю) внедиагональной связи в каждой из суперэлементных матриц жесткости, соответствующих блоку Л ^. Важной собенностью разработанной методики является то, что при использовании специального упорядочения неизвестных матрицы , полученные с помощью алгоритма, имеют блочно-диагональную структуру и могут быть обращены с арифметическими затратами, пропорциональными их размерности. Число обусловленности к ^Ау^Му? ^ при этом ограничено независимо как от параметров сетки, так и от анизотропности коэффициента эллиптического оператора:

«(4МГ1) <

1+\ 7 15

7 15

■5.31.

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

(к)

Рассмотрим предлагаемую методику построения блоков Мп подробнее.

Рассмотрим два последовательных уровня к — 1 и к (к — ко + 1,..., ¡) равномерного измельчения сетки, соответствующих триангуляциям 0к~ 1 и Г2ь Предположим, что процедура измельчения сетки такова, что каждый треугольник сетки содержит четыре конгруэнтных треугольника сетки П*. В дальнейшем мы будем называть совокупность этих четырех треугольников "суперэлементом" сетки Пь

Дизассемблируем билинейную форму Л(и,и), определенную выше, по отношению к суперэлементам и рассмотрим блоки Л^ суперэлементных матриц жесткости. Каждый блок А\^ соответствует "новым" узлам в рамках суперэлемента р (узлам из 14/^-1)- Очевидно, что глобальная матрица Лможет быть представлена посредством ассемблирования локальных

7 (к)

матриц АПр, а именно:

= Е

Р=1

где Л^ обозначает количество суперэлементов на уровне к, а матрицы Ц® соответствуют отображению глобальной нумерации узлов сетки на локальную нумерацию узлов в рамках суперэлементов. Без ограничения общности достаточно рассмотреть только случай произвольной формы суперэлемента и единичной матрицы коэффициентов [а,^]. В этом случае матрицы могут быть выписаны в явном виде:

-0Р 1 -Ор

№ -ЛН,р — -ар ар + /Зр + % II 1 -ъ

-0Р ~Ъ О!р+0р + Ъ . .-4 -ъ 1

где

?]р = ар + 0р + %, ар = 0р = , % =

Vр Щ Щ

а параметры ар, ¡Зр и % являются котангенсами соответствующих углов ар,

0Р и тр суперэлемента.

Предположим, что углы ар, ¡Зр и тр таковы, что |йр| > шах(|Д,|, |7,,|).

Определим в этом случае переобуславливатель (аппроксимацию) к

локальной матрице Аи'р как

1 ~ар 0 "

II с. 1 0

0 0 1

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

> тах(|йр|, |тр|) или \%\ > тах(|йр|, |Д,|),

аппроксимация М\\р строится аналогично, путем сохранения соответствующих внедиагональных элементов в Далее мы покажем, что матрицы

М$р и Л^р являются спектрально эквивалентными друг другу с константами спектральной эквивалентности, не зависящими от параметров сетки. Рассмотрим обобщенную задачу на собственные значения:

Собственные значения Ар определяются из условия

¿еь(лЦр-ЛрМЙ)= 0.

Поскольку матрицы Лп^, и М^^ имеют размерность 3, собственные числа Ар могут быть найдены в явном виде:

Ац» = 1,

= а Ш1К,р = 1 + \[/(а!р,Рр,7р),

-Чр = Атш,р = 1 -

где

1. 1Хр

Как может быть легко показано, для суперэлемента произвольной формы с углами ар, 0р и -ур такими, что \ар\ > тах(|/Зр|, |7Р|), функция /(йр,/Зр,7р) ограничена по модулю:

|/(йРЛ>7Р)1 <

из чего следует, что

тах,р ^ 1 \

7_ 15"

к (м«-1 Л\1) = < и 5.31.

У1Б

Следовательно,

1 + .

Очевидно, что такая же оценка справедлива и для случаев, когда |Д,| > шах(|ар|, |тр|) или \%\ > тах(|ар|, |Д>|),

« {и)

если только ненулевые связи в М11р будут выбраны соответствующим образом. Этот результат сформулирован в следующей теореме.

Теорема 2.1. Пусть тензор диффузии [а,^] —произвольная симметричная положительно определенная матрица, и пусть 01 — произвольно выбранная триангуляция. Тогда для всех матриц Л11р можно построить спек- /¡л — /1_\ Т

трально эквивалентные аппроксимации = Мп'р , содержащие только по два внедиагопальных элемента, такие, что справедлива следующая оценка:

г 7

0< 1-

1 +

М5-

Теорема 2.2. Представим матрицу Л/Ц как результат ассемблирования суперэлементных матриц

< = е

р=1

При использовании алгоритма построения локальных матриц Ми]р, основанного на сохранении единственной (максимальной по модулю) внедиаго-нальной связи в суперэлементной матрице жесткости Лп^, число обусловленности к Ли^ ограничено независимо от матрицы коэффициентов и параметров дискретизационной сетки. Более того, алгебраическая система с матрицей М$ может быть решена с арифметическими затратами, пропорциональными ее размерности.

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

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

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

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

Рассмотрим краевую задачу

Си = /, С = -У/х(х)У, ¿¿(х) >ц0>0, хеП

со стандартными граничными условиями, заданную в призме П, являющейся тензорным произведением П = Пг ® 0.ху отрезка числовой оси Г22 и многоугольника О.ту. Дискретизация данной задачи будет проводиться на иерархической последовательности конечно-элементных сеток

т М =%® Т«, 1 = 0 ,...,1,

где Тг — одномерная сетка, покрывающая £1г, а Т.М — последовательность вложенных триангуляции, покрывающих £1ху. Из определения Т^ следует, что она состоит из прямоугольных призм, ориентированных вдоль оси 2. В дальнейшем мы не будем налагать никаких ограничений на сетку Тг, но будем предполагать, что сетки Т^ не содержат треугольников, близких к вырожденным. Несмотря на указанное ограничение, рассматриваемый класс сеток допускает значительное варьирование шага дискретизации Л, что позволяет достигать хорошей аппроксимации задачи при использовании сеток со сравнительно небольшим количеством узлов.

Введем следующее естественное расщепление оператора задачи: С=Сжу + Сг, Сху = А^эдА, =-A^jL.

В дальнейшем мы будем предполагать, что функция /z(x) постоянна на каждой из призм, формирующих сетку Т'0', и может иметь разрывы на границе между ними. Мы предположим также, что оператор £ху дискрети-зирован посредством конформных кусочно-линейных конечных элементов, а оператор Сг — посредством конечных разностей. Рассмотрим последовательность нодальных матриц жесткости {А^ = + /Щ}, порожденных дискретизацией оператора С = Cz + Cxv на последовательности сеток {тЩ, а также последовательность матриц жесткости

{Д( о = Д(о + АЩ},

соответствующих иерархическому базису конечных элементов, заданных на последовательности сеток {Т^}. Матрицы А^ и связаны следующим образом:

где матрицы J^ содержат связи между узлами {Т^}, соответствующими разным уровням дискретизации. Заметим, что объем памяти, требуемой для хранения матриц J^ и JM , а также арифметические затраты, необходимые для их умножения на вектор, пропорциональны размерности матрицы Для решения линейных систем с матрицей AW предлагается использовать аддитивную версию алгебраического метода многоуровневых итераций. При этом для построения аппроксимаций М\х предлагается использовать тот же подход, что и в двумерном случае. Рассмотрим этот вопрос подробнее. Расщепление оператора L порождает аналогичное расщепление

3(0 _ 3« . 3(0 Аи — Aur,

где блоки Ay1 соответствуют оператору Сху, а блоки Л 1ц — оператору £z. Рассмотрим построение аппроксимаций

wtf = wil + й®

к блокам Ап.

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

Определим блок как

Очевидно, что матрицы Л^ = -<4пх, + и Й^ = + Й^ являются спектрально эквивалентными друг другу; при этом константы спектральной эквивалентности не зависят от параметров сетки Тг.

Благодаря структуре оператора С2, неизвестные в матрице связаны между собой только вдоль оси г?. Следовательно, возможно построение такого переупорядочения неизвестных, в котором матрица Йп' будет блочно-диагональной с блоками, являющимися трехдиагональными матрицами. Таким образом, решение систем с Й^ может проводиться с использованием быстрых прямых методов. Определим блок М^ следующим образом:

-1

где PVi(x) — чебышевский многочлен, степень которого выбирается из условия

«(Aiff 1ЛЙ) <К< 16-2= ю 1.15.

Поскольку матрицы А^ и спектрально эквивалентны друг другу, для любого наперед заданного К существует такое i'¡, что указанное условие удовлетворяется. Следовательно, при подходящем выборе V{ и коэффициентов Р^(х) рассматриваемый многосеточный метод имеет оптимальный порядок арифметической сложности.

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

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

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

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

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

5. Создан пакет программ, реализующих различные варианты алгебраических многосеточных методов, в том числе и на параллельных ЭВМ.

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Axelsson О., Georgiev К., Mellaard М., Neytcheva М., Padiy A., Scalable and optimal iterative solvers for linear and nonlinear problems. - preprint No 9613, University of Nijmegen, 1996.

2. Axelsson O., Padiy A., On the additive version of the algebraic multilevel iteration method for anisotropic elliptic problems. - preprint No 9711, University of Nijmegen, 1997, принято к печати в SIAM J. Sei. Comput.

3. Axelsson О., Padiy A., On a robust and scalable linear elasticity solver based on a saddle point formulation. - preprint No 9723, University of Nijmegen, 1997, принято к печати в Int. J. Num. Meth. Eng.

4. Padiy A., Multilevel iterative method for anisotropic elliptic problems // In: Proceedings of the Conference on Preconditioned Iterative Solution Methods for Large Scale Problems in Scientific Computations "PRISM-97", University of Nijmegen, May 27-29, 1997, pp.145-169.

5. Padiy A., On a robust multilevel method applied for solving large-scale linear elasticity problems. - preprint No 9806, University of Nijmegen, 1998, принято к печати в Comm. Num. Meth. Eng.

6. Ларин M. P., Падий А. В., Модельный анализ приближенного варианта метода окаймления // Труды Конференции Молодых Ученых, Новосибирск, 20-22 апреля, 1998, Издательство ИВМиМГ СО РАН, Новосибирск.

7. Padiy A., On a parallel multilevel solver for linear elasticity problems. -HIPERGEOS preprint, 1998.

8. Padiy A., Larin M., Model analysis of a subspace correction technique for anisotropic diffusion problems. - preprint No 9818, University of Nijmegen, 1998.

Текст научной работы диссертации и автореферата по математике, кандидата физико-математических наук, Падий, Александр Викторович, Москва

РОССИЙСКАЯ АКАДЕМИЯ НАУК ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ

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

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

АЛГЕБРАИЧЕСКИЕ МНОГОСЕТОЧНЫЕ МЕТОДЫ ДЛЯ ЗАДАЧ С МАЛЫМ ПАРАМЕТРОМ

Падий Александр Викторович 01.01.07 — Вычислительная математика

Научный руководитель: проф. Кузнецов Ю. А.

осква - 1998

Оглавление

Введение 3

1 Алгебраические методы многоуровневых итераций (АММИ) 9

1.1 Двухуровневый метод. Усиленное неравенство Коши-Буняковского..........9

1.2 Многоуровневый метод. Стабилизация числа обусловленности ..............14

2 Аддитивная версия АММИ для анизотропной задачи диффузии 16

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

2.2 Аддитивная версия АММИ с внутренними итерациями......................17

2.3 Обобщенный метод сопряженных градиентов..................................21

2.4 Алгоритм построения аппроксимаций к блокам A\i ..........................24

2.5 Численные эксперименты..........................................................32

3 Модельный анализ проекционного переобуславливателя 40

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

3.2 Модифицированный проекционный переобуславливатель......................42

3.3 Алгоритм построения блока V....................................................46

3.4 Численные эксперименты..........................................................47

4 Алгебраический многосеточный метод для трехмерных задач, дискре-тизированных на сетках, состоящих из призм 53

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

4.2 Оценка параметра а................................................................55

4.3 Алгоритм построения блоков ................................................57

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

5 Применение АММИ для решения задач теории упругости 61

5.1 Постановка задачи..................................................................61

5.2 Стандартная конечно-элементная дискретизация..............................64

5.3 Дискретизация на основе смешанной формулировки ..........................65

5.4 Численные экперименты ..........................................................70

Введение

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

Использованию итерационных методов в вычислительной практике посвящена обширная литература. Подробное изложение теории итерационных методов можно найти, например, в работах [14], [17], [21], [25]. Простейшими итерационными методами являются метод Ричардсона, метод Якоби, метод Гаусса-Зейделя, а также метод последовательной верхней релаксации. Для класса алгебраических систем, порождаемых конечно-элементной дискретизацией эллиптических уравнений второго порядка, скорость сходимости упомянутых выше итерационных методов зависит от параметра дискретизации к (размера шага дискретизационной сетки) и значительно ухудшается при стремлении И, к нулю. Скорость сходимости основных итерационных методов можно существенно улучшить, применив к ним процедуры ускорения, такие как чебышевское ускорение или ускорение по методу сопряжённых градиентов. У каждой из упомянутых процедур ускорения есть свои достоинства и недостатки. Реализация чебышевского итерационного метода более проста, при его использовании не требуется вычисления скалярных произведений, что наиболее ценно при проведении рассчетов на параллельных компьютерах. С другой стороны, существенным недостатком чебышевского итерационного метода является необходимость иметь достаточно хорошие оценки границ спектра матрицы, которые, как правило, недоступны. Поэтому метод сопряжённых градиентов является одним из наиболее привлекательных методов для решения систем линейных алгебраических уравнений с симметричными положительно определёнными (и полуопределёнными) матрицами.

И для чебышевского метода, и для метода сопряжённых градиентов скорость сходи-

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

||х - XfclU <2ак ||х - х0||а, ||х - х*||2 < 2 ||х - х0||2,

где х — точное решение, а = (л/к - 1)/(у/к+ 1), а к = Amax(A)/Amin(A) — спектральное число обусловленности матрицы системы. Аналогичная оценка справедлива также и для чебышевского итерационного метода:

||x-Xfc||2 < 2ak ||х — Хо||г-

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

Ах -Ъ

рассматривается система

М~1Ах = М~гЬ (либо, альтернативно, система АМ~ху — b, х = М-1 у)

с матрицей М, выбранной таким образом, что, с одной стороны, обращение M легко реализуемо и, с другой стороны, число обусловленности матрицы М~гА значительно меньше числа обусловленности исходной матрицы А. Для достижения скорости сходимости, не зависящей от параметров рассматриваемой алгебраической задачи (в том числе не зависящей от шага h дискретизационной сетки) необходимо построить спектрально эквивалентный переобуславливатель M такой, что к(М~1А) < Const. Для обеспечения приемлимых арифметических затрат при решении возникающих в современной вычислительной практике систем большой размерности (где количество неизвестных достигает сотен тысяч и даже миллионов) желательно также, чтобы стоимость умножения на переобуславливатель росла не более чем линейно с ростом размерности задачи.

При решении самосопряжённых эллиптических задач одним из возможных способов построения переобуславливателей, удовлетворяющих указанным выше требованиям, является использование многосеточных методов. Впервые многосеточный метод был предложен Р. П. Федоренко в работах [19] и [20], где была рассмотрена итерационная процедура решения уравнения Пуассона и был доказан оптимальный порядок ее вычислительной сложности. Исследование скорости сходимости многосеточного метода для разностных эллиптических операторов второго порядка с переменными коэффициентами было проведено Н. С. Бахваловым [2]. Для случая конечно-элементной аппроксимации эллиптического уравнения сходимость многосеточного метода была доказана Н. Г. Астраханцевым [1]. Среди первых работ, посвященных данной тематике, отметим также работы [13], [15], где В. И. Лебедевым был предложен KP-метод, близкий по своей структуре к многосеточным методам. В дальнейшем интерес к многосеточным методам начал стремительно нарастать. Упомянем следующие работы, внесшие существенный вклад в развитие теории многосеточных методов: W. Hackbusch [63], [64], [65]; A. Brandt [57], [58]; D. Braess [49], [50], [51]; R. Е. Bank, T. Dupont [43], [44], [45]; S. McCormick [77], [78], [79]; J.-F. Maître,

F. Musy [72], [73]; H. Yserentant [93], [94]; В. В. Шайдуров [22], [23], [24]; J. H. Bramble, J. E. Pasciak, J. Xu [55], [56]; O. Axelsson, P. Vassilevski [37], [38]; Ю. А. Кузнецов [66], [68], [69], [74]. Систематическое изложение теории многосеточных методов может быть найдено в монографиях [24], [52], [65], [88].

При исследовании скорости сходимости "классических" многосеточных методов существенно используется предположение о регулярности рассматриваемой дифференциальной задачи. Как было показано в работах [42] и [49], если дифференциальная задача имеет решение из Н2 для любой правой части из L^ то при использовании многосеточного метода типа F-цикла в качестве переобуславливателя итерационный процесс является сходящимся, при этом скорость сходимости не зависит от числа используемых уровней (от размерности решаемой алгебраической системы). Этот результат был в дальнейшем обобщён в [53], [60] на случай //^"-регулярных задач, 0 < а < 1. Независимость скорости сходимости многосеточного метода по типу W-цикла от числа уровней была доказана в работах [43] и [65] в предположении Нх +а-регулярности решения.

Одним из подходов к построению многосеточных методов, не опирающихся на какие-либо предположения о регулярности задачи, является построение переобуславливателей, основанных на двухуровневом разбиении сеточных функций. Алгоритм такого типа был рассмотрен в [26], [27], [72], а также в ряде других работ. На основе метода, предложенного в этих работах, Н. Yserentant [93] перешёл от двухуровневого разбиения к многоуровневому, т.е. к многосеточным переобуславливателям. В данном подходе (а также в его модификации, предложенной в [44]) используется метод конечных элементов с так называемыми иерархическими базисными функциями. Для метода конечных элементов с использованием классического "нодального" базиса многосеточный метод такого типа был предложен в работе [89]. Эти методы являются близкими к оптимальным для двумерных задач (в том смысле, что число обусловленности матрицы М~1А растёт как log Л,-1 или log2 h~l при h —> 0). Однако, они не являются достаточно эффективными при решении трёхмерных задач: число обусловленности растет как h~l.

Подход к стабилизации числа обусловленности данного класса многосеточных переобуславливателей был разработанн в работах [37], [38] (О. Axelsson и P. Vassilevski). Доказательство сходимости предложенного в этих работах алгебраического метода многоуровневых итераций (АММИ) опирается на использование алгебраического аналога так называемого усиленного неравенства Коши-Буняковского; стабилизация числа обусловленности достигается посредством рекурсивного использования в методе внутренних чебышевских итераций. Важной особенностью данного переобуславливателя является то, что он оптимален (или близок к оптимальному) по отношению как к скорости сходимости итерационного процесса, так и к вычислительным затратам на его реализацию; при чем это справедливо вне зависимости и от регулярности задачи, и от разрывности коэффициентов эллиптического оператора. Заметим также, что сходный многосеточный переобуславливатель, стабилизированный с помощью внутреннего чебышевского итерационного процесса, был независимо предложен Ю. А. Кузнецовым [9], [68], [69]. Интересной особенностью данного метода (MGDD переобуславливателя) является то, что его можно рассматривать и как многоуровневый метод разделения области [8], [70], и как многосеточный метод.

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

переобуславливателей такого типа является одной из главных целей проводимых нами исследований. Настоящая работа посвящена построению итерационных алгоритмов, скорость сходимости которых слабо зависит от "малого параметра" в задаче. Примером такого "малого параметра" может быть либо минимальный угол в используемой триангуляции, либо коэффициент анизотропии в тензоре диффузии, либо степень "вытяну-тости" области, в которой поставлена задача. При построении методов, устойчивых по отношению к "малому параметру" нами будет существенно использована техника, разработанная в [26], [37] и [38]. В связи с этим предлагаее алгоритмы могут быть отнесены к классу алгебраических методов многоуровневых итераций (АММИ), предложенных в указанных выше работах.

Настоящая диссертация состоит из Введения и пяти глав. В Главе 1 изложены основные подходы к постоению АММИ переобуславливателей. В частности, рассматриваются мультипликативный и аддитивный способы построения двухуровневого пе-реобуславливателя, а также подходы к стабилизации скорости сходимости метода в многоуровневом случае. В Главе 2 предлагается специальная модификация аддитивного АММИ, устойчивая по отношению к анизотропии. Одной из новых идей, используемых при построении метода, является особый выбор аппроксимаций В^ к блокам А^, возникающих при 2x2 блочном разбиении матриц жесткости А^ на каждом из уровней дискретизации:

А\ о I } "новые" базисные функции _

} "старые" базисные функции ' о, ■ ■ ■,

АЮ =

41 л12

А(к) А{к) Л21 л22

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

Как будет показано в дальнейшем, число обусловленности к при этом огра-

ничено независимо как от параметров сетки, так и от анизотропности коэффициента эллиптического оператора:

ав иМР Ч < -Ь=«5.31.

Необходимо подчеркнуть, что предлагаемый алгоритм позволяет конструировать матрицы ВЙ} автоматически, при этом не накладывается никаких ограничений на выбор триангуляции самого грубого уровня. Это является основным преимуществом разработанного подхода по сравнению с подходом, предложенным в [75]. Разработанный алгоритм может быть также использован в случае билинейной конечно-элементной дискретизации на декартовых сетках. И в этом случае, сохраняя только самые сильные связи в локальных суперэлементных матрицах жесткости, можно построить легко обращаемые спектрально эквивалентные аппроксимации к блокам . Другой новой идеей, используемой

при построении устойчивого по отношению к анизотропии многосеточного метода, является использование проекционного переобуславливателя при решении задач на самой грубой сетке. Этот вопрос будет подробно рассмотрен в главе 3, где будет предложен эвристический алгоритм решения существенно анизотропных задач диффузии. В основу предлагаемой методики положена техника окаймления [18], [33] со специальным выбором окаймляющих векторов. Заметим, что предлагаемая итерационная процедура по своей сути близка к так называемому "дефляционному" методу сопряженных градиентов, предложенному в [81], но имеет при этом более низкую арифметическую сложность и лучшие параллелизационные свойства. В Главе 4 разработан легко параллелизуемый переобуславливатель для трехмерных эллиптических задач второго порядка, дискрети-зированных на конечно-элементных сетках, состоящих из прямоугольных призм. Особенностью предлагаемого алгоритма является специальное расщепление оператора задачи, используя которое удается построить многосеточный итерационный процесс, скорость сходимости которого ограничена независимо от разрывов в коэффициенте задачи. Глава 5 посвящена применению разработанных многосеточных переобуславливателей к задачам теории упругости.

В диссертации получены следующие основные результаты:

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

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

• Построенные переобуславливатели успешно применены для решения задач теории упругости. Исследовано два типа дискрети