CC BY
140
ПАРАЛЛЕЛЬНЫЕ МЕТОДЫ И ТЕХНОЛОГИИ
___ __ __ _______К» 1
ДЕКОМПОЗИЦИИ ОБЛАСТЕЙ1
В.П. Ильин
Рассматриваются параллельные методы декомпозиции областей для решения трехмерных сеточных краевых задач, получаемых в результате конечно-элементных или конечно-объемных аппроксимаций. Данные проблемы являются «узким горлышком» среди различных этапов математического моделирования, поскольку современные требования к разрешающей способности сеточных алгоритмов приводят к необходимости решения систем линейных алгебраических уравнений с числом неизвестных в сотни миллионов и с очень плохой обусловленностью, что вызывает экстремальную ресурсоемкость расчетов. Описываются многопараметрические варианты алгоритмов с различной размерностью декомпозиции — одномерной, двумерной и трехмерной, — с пересечением или без пересечения подобластей, при использовании величин перехлеста как оптимизирующих параметров, а также с различными видами внутренних условий сопряжения на смежных границах (Дирихле, Неймана или третьего рода). Исследуются вариационные итерационные процессы крыловского типа в пространствах следов с разными предобуславливающими подходами: операторы Пуанкаре-Стеклова, блочный метод Чиммино, альтернирующий метод Шварца аддитивного типа, а также грубо-сеточная коррекция, являющаяся в определенном смысле упрощенным вариантом алгебраического многосеточного подхода. Проводится сравнительный анализ критериев эффективности распараллеливания на многопроцессорных вычислительных системах.
Ключевые слова: декомпозиция областей, трехмерные краевые задачи, сеточные аппроксимации, параллельные итерационные алгоритмы в пространствах Крылова, предобуславливающие операторы.
Введение
Разбиение сложной задачи на более простые подзадачи традиционно является распространенным подходом к построению эффективных численных методов решения многомерных задач математического моделирования, описываемых системами дифференциальных или/и интегральных уравнений, а также эквивалентными вариационными постановками. Здесь можно выделить такие структурные принципы, как расщепление по координатным осям (неявные методы переменных направлений), расщепление операторов (алгоритмы дробных шагов), расщепление по физическим процессам, а также декомпозиция расчетной области на подобласти [1-6]. Эти же подходы, из которых мы остановимся только на последнем — методе декомпозиции областей (МДО), — лежат в основе достижения масштабируемого распараллеливания на современных многопроцессорных — многоядерных вычислительных системах (МВС), в том числе гетерогенной архитектуры, когда каждый вычислительный узел содержит один или несколько графических процессорных элементов (ГПЭ) общего назначения, обладающих высокой скоростью выполнения арифметических операций.
С точки зрения эффективного распараллеливания, необходимо рассмотреть типовую структуру алгоритмов, реализуемых в крупномасштабных вычислительных экспериментах посредством различных вложенных циклов. К самому внутреннему уровню следует отнести решение систем линейных алгебраических уравнений (СЛАУ)
1 Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2012».
высокого порядка (десятки или сотни миллионов), которое осуществляется с помо-тттью предобусловленных крыловских итерационных процессов. Эти «линейные» итерации могут иметь многоуровневый характер, когда СЛАУ имеют крупноблочную структуру, что имеет место при сеточной аппроксимации систем исходных функциональных уравнений с неизвестными векторными функциями. Следующий уровень — проведение нелинейных итераций, если коэффициенты и правые части алгебраической системы зависят от искомого решения. Для нестационарных задач предыдущие операции выполняются в цикле по временным шагам. Все указанные выше действия составляют реализацию прямых задач, а в случае решения обратных задач с применением оптимизационных подходов осуществляется внешний цикл с последовательным решением прямых задач и корректировкой целевого функционала на каждой итерации. «Узким горлышком» являются как раз самые простые — линейные — задачи, представляющие наибольший интерес с точки зрения распараллеливания, поскольку требуемый ими объем вычислений нелинейным образом растет при увеличении размерности, или числа степеней свободы.
Итоговая производительность программной реализации численных методов, конечно, значительно зависит от технологических аспектов, возникающих на различных стадиях наукоемких расчетов: геометрического и функционального моделирования, включающего интеллектуальные пользовательские интерфейсы, генерации адаптивных сеток (дискретизации расчетной области), аппроксимации исходной задачи методами конечных элементов, конечных объемов (МКЭ, МКО) и т.д., решения возникающих алгебраических систем сверхвысоких порядков, реализации оптимизационных процедур (в обратных задачах), а также постпроцессинга и визуализации результатов, см. [7]. Параллельные алгоритмы декомпозиции и формирование соответствующих структур данных, распределенных по вычислительным узлам МВС, должны рассматриваться комплексно на всех указанных этапах математического моделирования.
Однако главной целью данной работы является обзор и сравнительный анализ итерационных алгебраических методов декомпозиции областей, базирующихся, во-первых, на эффективном предобуславливании исходных сеточных СЛАУ и, во-вторых, на построении крыловских алгоритмов в пространствах следов, т.е. для функций и операторов типа Пуанкаре-Стеклова, определенных на внутренних границах смежных подобластей.
Помимо приводимых нами в списке литературы основных монографий [3-6] и некоторых статей [8-17] (отметим здесь значительный вклад новосибирских и московских математиков), по данной теме огромное количество материала имеется на интернетовском сайте [18], включая труды уже состоявшихся 20 специальных конференций.
1. Математические и технологические вопросы декомпозиции областей
Представим общую идею МДО следующим образом. Пусть в расчетной трехмерной области П с границей Г требуется решить краевую задачу
Ьи = f (г), г е П; 1и|г = д, (1)
где линейные дифференциальные операторы Ь и I обеспечивают существование единственного решения „(г) в некотором функциональном пространстве. Разобъем область П на Р пересекающихся или не пересекающихся подобластей Пд и введем следующие обозначения:
р
п=ипд, П = ^Г, Пд = П^Гд,
9=1 _ (2)
Гд = и Гд,д' > Гд,д' = Гд П Пд' > = 5>
д'&Шд
где шд есть множество номеров подобластей, соседних к Пд, Гд — вся граница q-й
подобласти, а Гд0 — ее внешняя часть, т.е. П0 обозначает «внешнюю» подобласть —
дополнение к П в Вместо исходной задачи (1) рассматриваем Р краевых подзадач в подобластях Пд:
Ьид (г ) = fд , Г е Пд ,
1д,д' (ид )|Гд,д/ = дд,д/ = 1д',д (ид/) I Гд/д ’ П9 |Гд,д' = ид/ |Гд>д/, (3)
q/ е ^д, /д.0ид|Гд,о = д, q =1,...,р,
где 1д,д/ и ддд/ при q/ = 0 определяют некоторые «внутренние» граничные условия такие, которые не «портили» бы исходное решение и, т.е. ид(г) = и(г) при г е Пд. Например, в достаточно общем случае эти условия имеют вид
, а ди9| , а дид/ 1
адид + вд о 1 г , дд,д/ — ад/ ид/ + вд/ о т , ,
дпд дпд/ 'V(4)
К| + |вд| > 0, ад ■ вд > 0,
причем для вд = вд/ = 0 они соответствуют условию 1-го рода (Дирихле), при ад = ад/ = 0 — условию 2-го рода (Неймана), а в остальных случаях — условию 3-го рода (Робина). Естественным методом решения системы краевых задач (3) является организация простых итераций, т.е. блочного метода Якоби вида
ЬиП = Л, 1д,д/иП|ГдУ = ^д^/-1!^ , (5)
где правые части граничных условий для каждой подобласти определяются по значениям решения с прошлой итерации в смежных подобластях.
Далее в качестве иллюстрации исходной задачи (1) мы будем рассматривать аппроксимацию уравнения Пуассона в кубической области с граничными условиями Дирихле на кубической сетке с общим числом узлов М3:
—Ди = Л, и|г = д, П = [0 х 1]3; П = {г, ^ к},
(Аи^)г)^;Й = 6и^ — и^-1)^;й — и^-1,А;-
? _ .
иі+1,І,к аі,і,к— 1 аі,і,к+1 ^і,І,к.
і,л* = і,...,м, /? = {/,}, и? = каіек'
Рис. 1. Примеры структурированных Ш-, 2Р- и ЗР-декомпозиций области
Отметим, что с точек зрения как скорости сходимости итераций, так и технологии распараллеливания, существенное значение имеет размерность декомпозиции, варианты которой представлены на рис. 1. Для простой расчетной области в форме параллелепипеда размерность очевидным образом определяется как количество типов координатных плоскостей, используемых при построении подобластей. Например, при 1-Р декомпозиции область разбивается на Р слоев с помощью параллельных плоскостей.
2. Итерационные методы в пространстве следов
Рассмотрим основные алгоритмические вопросы МДО на примере одномерной декомпозиции с пересечением подобластей описанной в (6) сеточной краевой задачи, см. рис. 2.
—I_______________________|_
1—1—1______I____. 1.
1—1—1______I____І І I
_|________________________|_
к] к, кГ к,
Рис. 2. Пример одномерной декомпозиции с пересечениями подобластей
Через Qh обозначаем сеточную подобласть, у которой первая и последняя координатные плоскости г = имеют номера к1 и кN соответственно. Все подобласти, как и их пересечения Ад, считаем одинаковыми. Введем также следующие обозначения:
N = ^гт(П^), т = ^гт(А^), д = 1,..., Р,
N = к^ - к? + 1, т = к^ - к?+1 + 1, М = PN - (Р - 1)т, (7)
Здесь N и т означают число г-плоскостей в Qq и Ад, М есть «размерность по г» всей сеточной области П^,а Р — общее количество подобластей.
Вводя подвекторы и? = (и*л,...,и*л )Т, Пц = {щ^ц; г,.?, = 1,...,М} е Пм2, а также аналогичные подвекторы правых частей /д порядка М2N, системы уравнений для подобластей можно записать в следующем блочно-трехдиагональном виде:
Ад,д- 1 ид- 1 + ид Ад,д+1 ид+ 1 /д, ? 1 •••, Р>
(8)
Л — Л _пл Л — ЛТ с им 2ы,м2 N
А1,0 = АР,Р+1 = ° лд,д, Ад,д±1 = Лд±1>д е .
Итерационный процесс (5) при использовании значений и^.-1 с предыдущего шага из соседних подобластей имеет вид
А ип = Г1-1 = / + Г1-1 + Г1-1
Ад,д ид = / — / ,
^ (9)
т-1 = а ,ип-1 Гп-1 = А ,ип-1
/ д = Ад,д-1ид-1, / д = Ад,д+1 ид+1 ,
который в по-компонентной форме записывается следующим образом:
(С - 0/- и%-1 = /к? + г#-1,
^пг! = ип-_11 - яип-1, к = кд,
(Aq,q un)fc — <
(С - вІ)unN - unN+1 = ffcN + Ч+l1, wn+l1 = un-+l - в<-1, k = kN,
—Un-1 + Cun — Un+1 = ffc, k k = kf + 1, •••, kN — 1 •
Здесь I — единичная, а С — пятидиагональная матрицы порядка M2, т.е.
(С — вІ )un)fc )i,j = {(б — — un-l,j,fc — un+l,j,fc — Unj-l,fc — Un!j+l,fc
Un-1 Є ftq_
С ЄКМ2 , в Є [0, 1], Un-1 Є fiq-1,
а в — итерационный параметр, регулирующий тип итерируемого граничного условия на смежных границах подобластей: значения в = 0, в = 1 и 0 <в< 1 соответствуют условиям Дирихле, Неймана и Робина.
Вводя далее подвекторы и размерности M2, соответствующие лежащим в сеточным слоям с номерами 1 и &1+1, итерационные соотношения (9) приведем к редуцированному виду
v
"n - Bq,q- lW^n-1 + Bq,q+lvqn^1 + gq, q = 2, •••, P,
q _ q,q- 1 q- 1 1 -^q
wn = Bq,q-1 Wq-1 + ^q.q+l^+l + gq, q =1,•••,P — 1, (l0)
B1,0 = BP,P +1 = °
где используются следующие обозначения (см. подробнее [9]):
vq Cq,q—1Uq, wq Cq,q+1uq, Aq,q±1 Qq,q±1 C^qil,
I-1 f
q,q fq,
gq Cq.q+l А-,о/ї,
Вд,д±1 Сд,д—1 Вд,д±1 Сд,д+1 А—д^д±1 •
Система уравнений (10) получена из (9) исключением «внутренних» неизвестных для подобластей и содержит только подвекторы, соответствующие смежным границам. Объединяя их в «большие» векторы следов 5 = (^1, г>2,^р-1, ^р)т, д = (^,$2, ...,&-1 ,др)Т, размерности 2М2(Р — 1), итерации (10) перепишем в сжатом виде
= Тзга-1 + д, п = 1, 2,..., (11)
где явный вид матрицы Т нам не потребуется. Отметим только, что умножение на нее включает решение алгебраических подсистем для всех подобластей. Очевидно,
что если итерационный процесс (11) в пространстве следов сходится, т.е. Зп ^ 5, то предельный вектор удовлетворяет системе уравнений
А5 — (/ - Т)5 = 2, (12)
в которой А есть некоторая предобусловленная (по отношению к А из (6)) матрица.
Эффективное численное решение полученной предобусловленной СЛАУ реализуется с помощью какого-либо из итерационных алгоритмов в подпространствах Крылова, см. [19] и цитируемую там литературу. Например, если А есть симметрично положительно определенная матрица, то целесообразно применять методы сопряженных градиентов или сопряженных невязок, которые при V = 0,1 соответственно описываются следующими единообразными формулами:
г
0 - д — Аз0 = 51 — 50, 51 = Т50 + д, р0 = г0, п =1, 2,
5п+1 = 5п + Оп'У, = р^ )М1°, Рп} = (А^ гп, гп),
^) = (Арп, А(^рп), ) = рпЙ-1/рп1'),
гп+1 = гп - )^4рп, рп+1 = гп+1 + в^)Рп.
В более общем случае, когда А несимметрична, алгоритмы строятся на основе А^ ортогонализации Арнольди (V = 0 или V =1):
ип = и0 + у^1 + ... + Уп^п, К,А"^) = ^Т^п, ^ = (^п,А"^п),
п+1 = Vй , V1 = г0 = f — Аи0,
= (А8^,^)
л<> ) = (А^п,А^ к =1 П + 1 V -, = (?)1 ^™+1)
(Л8^ к 1,...,п + 1, ^га+1 (^ ,...,^ )
я» = {ЛІГІ}
я„
еП
є 7г»+1,», я» є 7г»
При этом коэффициенты у., к = 1,..., п, разложения итерационного приближения
ип по базисным векторам V1,..., V”- находятся из условий ортогональности или мини-
к «->
мальности векторов невязок г , что приводит к методам полной ортогональности или
обобщенных минимальных невязок (РОМ, А-ЕОМ, GMR.ES, A-GMR.ES).
3. Итерационные процессы на основе операторов Пуанкаре—Стеклова
В данном пункте мы рассмотрим частный случай декомпозиции на две подобласти без пересечения (Р = 2, А = 0), в которых решаются уравнения Пуассона
-Аид = /д, и|г = 2,9 = 1, 2, (13)
совместно определяющие в расчетной области П = П1 и П2 решение задачи Дирихле, обладающее свойствами непрерывности на общей границе подобластей Г1>2:
-Аи = /, и = и^ и2; Г 1,2 : и1 = и2, -. (14)
Задавая на внутренней границе произвольное начальное приближение и вычисляя решения соответствующих подзадач
u0 = и2|Г1,2 = ur 1,2, -Au0 = Л, uV = g (15)
для функций ошибки vq = uq — «0, q = 1, 2, мы получаем следующие уравнения:
-Д^9 = 0, v9|г = о, v9|ri,2 = иГ 1,2 — «Г1,2 = vr 1,2 • (16)
Отметим, что в силу условий непрерывности исходного решения, имеет место равенство
dvi + = = _/ диЦ + А (17)
5ni dn2 ^ V dn dn2/ •
Определяя теперь на Г12 операторы Пуанкаре-Стеклова Sq
dv
= S— 1vr 1,2 = Sg“1(ur1,2 - «Г12), q = 1, 2, (18)
мы приходим к уравнению
,4uru = (S-1 + S-1)uru = Ф = (S-1 + S-1)u?1J - (|П1 + M), (19)
которое фактически представляет новую формулировку метода декомпозиции областей. Очевидно, что умножение на оператор S-1 включает решение краевой подзадачи в .
Оказывается, что уравнение (19) может быть эффективно предобусловлено следующим образом:
Au = BAu = B^, B = S1 + S2. (20)
Предобусловленный оператор A обладает при этом замечательными свойствами: A = (S1 + S2 )(S-1 + S-1) = I + S1S-1 + I + S2S-1 =
= A1 + A2, A1 = I + S1S-1, A2 = I + S2S-1, A1A2 = A1 + A2 = A2A1,
т.е. он представим в виде суммы перестановочных операторов A1, A2. Это позволяет, в частности, применять для решения (20) сверхбыстрые неявные методы переменных направлений [1, 2, 19]. Отметим также, что в работе [9] на основе операторов Пуанкаре-Стеклова построен оптимальный (сходящийся за 2 итерации) метод Шварца при декомпозиции области на две непересекающиеся подобласти.
4. Блочный метод Чиммино
Фактически проективный подход к алгебраической декомпозиции реализуется в блочном методе Чиммино, см. [20] и цитируемую там литературу. Разбивая вектор правой части f = (Д,..., fp)т на подвекторы, мы можем записать СЛАУ в блочном виде
A1 ' f 1 '
Au — u =
Ap fp
— f, A gRn,n; u,f gR
(2l)
где Ад суть блочные строки матрицы, которые для простоты считаем содержащими одинаковое число М строк. Записывая соответствующие подсистемы в виде
Ади = /, / е^м, Ад , к = 1,..., Р, N = РМ, (22)
и вводя вспомогательные подвекторы
< = А+гД, гг = /к - Аки", А+ = АТ(АкАТ)-1, (23)
мы можем определить итерационный процесс
р
и"+1 = и" + ш^г;", п = 0,1,..., (24)
к=1
где ш есть некоторый итерационный параметр, а А+ — псевдообратная матрица, корректно определяемая, если Ак есть матрица полного ранга.
Легко видеть, что на каждом п-м шаге вычисляются («одновременно» для всех к) в некотором смысле приближенные решения в соответствующих подобластях.
Стационарный итерационный процесс (24), очевидно, может быть ускорен с помощью какого-то из крыловских методов. Рассмотренный алгоритм можно усилить, если в (23) А+ заменить на обобщенную псевдообратную матрицу
А+С = С-1АТ(АдС-1 АТ)-1, Ск е Ям'м, (25)
где Сд — некоторая симметричная положительно определенная матрица. Умножение на участвующую в (25) обратную матрицу может быть реализовано с помощью решения вспомогательной системы с седловой точкой
' С АТ! ' <' 0
Ак 0 . < . Г.П 'к
Сдо" = -АТV", V" = А+Сг".
Понятно, что существующий большой произвол в выборе матриц Сд дает потенциально значительные возможности в ускорении итераций.
5. Аддитивный метод Шварца с грубо-сеточной коррекцией
Если расчетную облалсть представить в форме П = Пд У Пя, где Пя — дополнение к подобласти Пд в П, а также ввести обозначения
и
ид
ид = Яд и,
(27)
где Яд = [0/0] — матрица продолжения, а ЯТ : и ^ ид — матрица сужения вектора, то итерационный метод декомпозиции записывается в виде
рр
ига+1 = ип
+ 5] Вд(/ - Аип) = ип + Вгга, В = ^ Вд,
д=1 д=1
который называется аддитивным алгоритмом Шварца. Здесь предобуславливающая матрица В состоит из слагаемых
Вд = ЯТ (Яд АЯТ )_1Яд, (29)
каждое из которых является симметричным при А = АТ и положительно определенным, если таким же свойством обладает А, а матрица Яд имеет полный ранг.
Матричное произведение Рд = ВдА обладает свойством Р2 = Рд, то есть является ортогональным проектором в смысле А-скалярного произведения
(Рди, 1>)а = иТР^А-у = иТАВдА-у = (и, Рд^)А.
Сходимость итерационного процесса (28) можно ускорить за счет какого-либо «улучшения» предобуславливателя В. Мы рассмотрим метод грубосеточной коррекции [4-6], который заключается в следующем.
Запишем исходную алгебраическую систему в виде
Арир = /р, ир, /р е ^, Ар е ^’^, (30)
предполагая, что она получена из аппроксимации некоторой краевой задачи на густой сетке Пр и имеет достаточно большую размерность Nр. Определим теперь подвектор
ис = Ясир, Яс , Nc < Nр,
который будем считать соответствующим некоторой редкой (грубой) сетке Пс, вложенной в Пр, а также дополнительный предобуславливающий оператор
Вс = ЯТА-1Яс, Ас = ЯсАЯТ е ^Мс’Мс.
Тогда для решения СЛАУ можно определить аддитивный метод Шварца с грубосеточной коррекцией следующего вида:
ир+1 = ир + (Вс + Вр )(/р — Ар ир), (31)
где Вр = ^р=1 Вд — новое обозначение предобуславливателя В из (28).
В частности, если матрица ЯТ соответствует оператору продолжения интерполяционного типа с сетки Пс на Пр, то получаемый алгоритм называется методом галер-кинского вида. Естественно, что итерационный процесс (31) может быть ускорен с помощью крыловских подходов. Введение дополнительного предобуславливателя Вс реализует на каждой итерации взаимосвязи между дальними подобластями и устанавливает некоторый «мостик» между декомпозицией областей и многосеточными методами.
В заключение данного пункта отметим еще такой независимый подход к ускорению итераций, как метод дефляции, см. [17] и приведенный там обзор. В применении к СЛАУ (30) его можно описать с помощью матрицы некоторого проективного пространства Я^ е с полным рангом и заданным к < Nр. Используя Я^, можно
определить проектор
Рй = 1 — арв, в = я^(яТаря^_1яТ
и решать вместо (30) предобусловленную им систему Р^Арир = Р/р. При определенных выборах Я^ (например, использование приближенных собственных векторов Ар) метод дефляции дает значительное ускорение итераций.
6. Вопросы распараллеливания методов декомпозиции
Основные критерии распараллеливания алгоритмов — это коэффициенты ускорения и эффективности использования процессоров Ер, которые выражаются через время выполнения арифметических операций на Р процессорах Тр и длительности межпроцессорных коммуникаций Тр:
= Т1/ТР, Ер = /Р,
(32)
ТР = Тр + Тр ~ + Дс(т0 + тсЮ-
Здесь V, означает общее количество арифметических действий, та - некоторое усредненное время выполнения одной операции, N — число обменов одного процессора, V — количество передаваемых чисел за одну коммуникацию, тс — время передачи одного числа, а т0 - длительность задержки при каждом обмене, причем можно считать та ^ тс ^ т0. При этом вынужденно используется грубая модель вычислительного процесса, которая не учитывает многих факторов современных МВС: гетерогенности архитектуры и неоднородности иерархической памяти, использования конвейеризации и совмещения во времени арифметических действий с обменами, топологии физической реализации межпроцессорных соединений, особенностей функционирования многоядерных вычислительных узлов и т.д. В отсутствие реальной возможности реального имитационного моделирования компьютерной системы теоретические оценки ее производительности могут иметь только качественный характер, и практически единственным средством исследования эффективности распараллеливания алгоритмов решения определенного класса задач является численный эксперимент.
Проведем качественный сравнительный анализ эффективности распараллеливания ё-Р декомпозиции трехмерной области для описанной в (6) модельной сеточной задачи для различных размерностей й = 1, 2, 3, изображенных на рис. 1. Пусть во всех рассматриваемых случаях область разбивается на Р одинаковых сеточных подобластей с числом узлов М3/Р. Количество же граничных узлов каждой подобласти равно
л'!-1’ = 2М2( 1 + 1). л'!-2’ = 2М2( ^ + )• N(3) = ^3^
для й = 1, 2, 3 соответственно. Поскольку при Р ^ 1 имеем N(1) ^ N(2) ^ N(3), а объем данных, передаваемых от одного процессора ко всем остальным, пропорционален числу поверхностных сеточных узлов, то очевидным образом получаем, что относительный вклад коммуникационных временных потерь уменьшается с увеличением размерности й.
Большую роль играет также величина, которую можно назвать псевдодиаметром ж графа макросети, образуемой совокупностью подобластей. Расстоянием между двумя вершинами графа называется минимальное число ребер, по которым можно пройти из одной вершины в другую. А псевдодиаметр — это максимальное расстояние между какими-либо парами вершин. Если кубическую сеточную область из (6) разбить на Р подобластей с помощью ё-Р декомпозиции, то псевдодиаметры будут равны ж = тДР 1/й, т.е. длине диагонали квадрата или куба при й = 2, 3 соответственно.
Очевидно, что при реализации блочного метода Якоби по подобластям информация от одной подобласти дойдет за одну итерацию только до соседних подобластей и дойдет до самой дальней подобласти за ж итераций в худшем случае. Можно сделать естественное предположение (теоремы такой пока нет), что число итераций по подобластям, необходимое для достижения требуемой точности е ^ 1, пропорционально ж7|/ne|, y > 0, т.е. при P > 8 (минимальное число подобластей в трехмерной декомпозиции), т.е. будет убывать с ростом размерности d.
Таким образом, с двух различных теоретических точек зрения, 3-D декомпозиция является самой предпочтительной. Однако на практическое соотношение производительности программных кодов, реализующих различные подходы, естественно, могут повлиять многочисленные технологические факторы. Для рассмотренной в п.3 одномерной декомпозиции на одной итерации времена выполнения арифметических действий и обменов для каждого процессора равны (C1 — не зависящая от М и N константа)
T“ = C'iM2N7+1ra,
= 2(то + тсМ2)
и при этом критерии эффективности распараллеливания данного алгоритма при N ^ 1 оцениваются величинами
Sp = Tp ■ P/(T£ + Tp) « P, Ep « 1,
т.е. обеспечивают примерно линейное ускорение с ростом P.
Как уже отмечалось, фактическая производительность программной реализации может сильно отличаться от теоретических оценок. В частности, для распараллеливания существенное значение имеет такая нетривиальная в общем случае сложных неструктурированных сеток задача, как построение сбалансированных по числу узлов подобластей, которую логично осуществлять на этапе дискретизации расчетной области. А вопросы оптимизации величин пересечений и типов внутренних граничных условий требуют дополнительных и аналитических, и экспериментальных исследований.
Работа поддержана грантом РФФИ №11-01-00205, а также грантами Президиума РАН №2.5 и ОМН РАН № 1.3-4-
Литература
1. Марчук, Г.И. Методы вычислительной математики / Г.И. Марчук. - М.: Наука, 1980.
2. Ильин, В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений / В.П. Ильин. - Новосибирск, изд. ИВМ СО РАН, 2000.
3. Лебедев, В.И. Операторы Пуанкаре - Стеклова и их приложения в анализе /
В.И. Лебедев, В.И. Агошков, - М.: Отдел вычислительной математики АН СССР, изд. ВИНИТИ, 1983.
4. Quarteroni, A. Domain Decomposition Methods for Partial Differential Equations /
A. Quarteroni, A. Valli - Clarendon Press, Oxford, 1999.
5. Smith, B.F. Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations / B.F. Smith, P.E. Bjorstad, W.D. Gropp - Cambridge University Press, 2004.
6. Toselli, A. Domain Decomposition Methods. Algorithms and Theory / A. Toselli,
O. Widlund - Springer, Berlin, 2005.
Т. Ильин, В.П. Параллельные процессы на этапах петафлопного моделирования / В.П. Ильин / / Вычислительные методы и программирование. - 2011.-Т. 12, № 1.
- С. 93-99.
8. Nataf, F. Optimized Schwarz Methods. // Lecture Notes in Computer Science and Engineering. - Springer-Verlag, Berlin, 2009. - P. 233-240.
9. Ильин, В.П. Параллельные методы декомпозиции в пространствах следов /
B.П. Ильин, Д.В. Кныш // Вычислительные методы и программирование. - Изд. МГУ, 2011. - Т. 12, № 1. - С. 100-109.
10. Смелов, В.В. Принцип итерирования по подобластям в задачах с эллиптическим уравнением. / В.В. Смелов В.В., Т.Б. Журавлева. - М.: Изд. ВИНИТИ, 1981. -(Препринт / ОВМ РАН; № 14).
11. Сандер, С.А. Модификация алгоритма Шварца для решения сеточных краевых задач в областях, составленных из прямоугольников и параллелепипедов. /
C.А. Сандер. - Новосибирск, 1981. - (Препринт / Изд. ВЦ СО АН СССР; № 83).
12. Мацокин, А.М. Применение окаймления при решении систем сеточных уравнений / А.М. Мацокин, С.В. Непомнящих // Вычислительные алгоритмы в задачах матемаической физики - Новосибирск: Изд. ВЦ СО АН СССР, 1983. - С. 99-109.
13. Лебедев, В.И. Вариационные алгоритмы метода разделения области / В.И. Лебедев, В.И. Агошков - М., 1983. - (Препринт / ОВМ РАН, № 54).
14. Непомнящих, С.В. О применении метода окаймления к смешанной краевой задаче для эллиптических уравнений и осеточных нормах в W21/2(S). / С.В. Непомнящих.
- Новосибирск, 1984. - (Препринт / Изд. ВЦ СОАН СССР, № 10б).
15. Кузнецов, Ю.А. Новые алгоритмы приближенной реализации неявных разностных схем / Ю.А. Кузнецов - М., 198Т. - (Препринт / ОВМ АН СССР, № 142).
16. Свешников, В.М. Построение прямых и итерационных методов декомпозиции /
B.М. Свешников // Сиб. журн. индустр. математики. - 2009. - Т. 12, № 3(39). -
C. 99-109.
17. Tang, J.M. Comparison of Two-level Preconditioners Derived from Deflation, Domain Decomposition and Multigrid Methods / J.M. Tang, R. Nabben, C. Vuik, Y.A. Erlangga // J. Sci. Comput. - 2009. - V. 39. - P. 340-370.
18. Domain Decomposition Methods.
URL: http://ddm.org (дата обращения: 14.03.2012)
19. Ильин, В.П. Методы и технологии конечных элементов / В.П. Ильин - Новосибирск, изд. ИВМиМГ СО РАН, 2007.
20. Ильин, В.П. Об итерационном методе Качмажа и его обобщениях. // Сиб. журн. индустр. математики. - 200б. - Т. 9, № 3. - С. 39-49.
Ильин Валерий Павлович, доктор физико-математических наук, профессор, главный научный сотрудник ИВМиМГ СО РАН, профессор кафедры вычислительной математики НГУ, ilin@sscc.ru.
PARALLEL METHODS AND TECHNOLOGIES OF DOMAIN DECOMPOSITION
V.P. Il’in, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation)
Parallel domain decomposition methods for solving 3-D grid boundary value problems, which are obtained by finite-element or finite-volume approximations are considered. These problems present the bottle neck between different stages of mathematical modelling, because the modern requirements to accuracy of grid algorithms provide the necessity of solving the systems of linear algebraic equations with the hundred millions of degrees of freedom and with super-high condition numbers which demand the extremal computing resourses. Multi-parameter versions of algorithms with various domain decomposition dimensions — one-dimensional, two-dimensional and three-dimensional, — with or without overlapping of subdomains and with different kinds of internal conjecture conditions on the adjacent boundaries (Dirichlet, Neuman and Robin). The iterative Krylov processes in the trace spaces are investigated for the different preconditioning approaches: Poincare - Steklov operators, block Cimmino method, alternating Schwartz algorithm of additive type, as well as coarse grid correction which is, in a sense, the simplified version of algebraic multigrid method. The comparative analysis of the criteria of parallelezation for the multiprocessor computer systems is made.
Ключевые слова: domain decomposition, tridimensional boundary value problems, grid approximations, parallel iterative algorithms in Krylov spaces, preconditioning operators.
References
1. Marchuk G.I. Metody vychislitelnoj matematiki [Numerical Analysis Methods] Moscow, Nauka, 1980.
2. Il’in V.P. Metody konechnyh raznostej i konechnyh ob”emov dlja jellipticheskih uravnenij [Finite Difference and Finite-volume Methods for Elliptic Partial Difference Equations]. Novosibirsk, Institute of Computational Modelling SB RAS, 2000.
3. Levedev V.I., Agoshkov V.I. Operatory Puankare - Steklova i ih prilozhenija v analize [Poincare-Steklov Operators and their Application in Analysis] Moscow, Computational Mathematics Department of USSRAS, VINITI, 1983.
4. Quarteroni A. Valli A. Domain Decomposition Methods for Partial Differential Equations. Clarendon Press, Oxford, 1999.
5. Smith B.F., Bjorstad P.E., Gropp W.D. Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations. Cambridge University Press, 2004.
6. Toselli A., Widlund O. Domain Decomposition Methods. Algorithms and Theory. Springer, Berlin, 2005.
7. Il’in V.P. Parallel’nye processy na jetapah petaflopnogo modelirovanija [Parallel Processes in Petaflopic Modeling]. Vychislitel’nye metody i programmirovanie [Numerical Methods and Programming], 2011. V. 12, No 1. P. 93-99.
8. Nataf F. Optimized Schwarz Methods. Lecture Notes in Computer Science and Engineering. Springer-Verlag, Berlin, 2009. P. 233-240.
9. Il’in V.P., Knysh D.V. Parallel’nye metody dekompozicii v prostranstvah sledov [Parallel Methods of Decomposition in Trace Spaces] Vychislitel’nye metody i programmirovanie [Numerical Methods and Programming]. MSU Publ., 2011. V. 12, No 1. P. 100-109.
10. Smelov V.V., Zhuravleva T.B. Princip iterirovanija po podoblastjam v zadachah s jellipticheskim uravneniem [Subdomain Iteration Principle in Partial Differential Equation Problems]. Moscow, VINITI, 1981.
11. Sander S.A. Modifikacija algoritma Shvarca dlja reshenija setochnyh kraevyh zadach v oblastjah, sostavlennyh iz prjamougol’nikov i parallelepipedov [Schwartz Algorithm Modification for Solving Grid Boundary Value Problems in Areas of Rectangles and Parallelepipeds. Novosibirsk, Preprint No 83, CC SB USSRAS, 1981.
12. Matsokin A.M., Nepomnyaschikh S.V. Применение окаймления при решении систем сеточных уравнений [Using the Bordering Method for Solving Systems of Mesh Equations]. Vychislitel’nye algoritmy v zadachah matematicheskoj fiziki [Numerical Methods in Mathematical Physics Problems]. Novosibirsk, CC SB USSRAS, 1981. P. 99-109.
13. Lebedev V.I., Agoshkov V.I. Variacionnye algoritmy metoda razdelenija oblasti [The Variational Algorithms of Domain Decomposition Method]. Moscow, Preprint No 54, Department of Numerical Mathematics of RAS, 1983.
14. Nepomnyaschikh S.V. O primenenii metoda okajmlenija k smeshannoj kraevoj zadache dlja jellipticheskih uravnenij i osetochnyh normah v W21/2(S)] [On the Application of the Bordering Method to the Mixed Boundary Value Problem for Elliptic Equations and on Mesh Norms in W2 (S)]. Novosibirsk, Preprint No 106, CC SB USSRAS, 1984.
15. Kuznetsov Yu.A. Novye algoritmy priblizhennoj realizacii nejavnyh raznostnyh shem [New Algorithms for Approximate Implementation of Implicit Difference Schemes] Moscow, Preprint No 142, Department of Numerical Mathematics of RAS, 1987.
16. Свешников В.М. Построение прямых и итерационных методов декомпозиции [Construction of Direct and Iterative Decomposition Methods]. Sib. Zh. Ind. Mat. [J. Appl. Industr. Math.], 2009. V. 12, No 3(39). P. 99-109.
17. Tang J.M., Nabben R., Vuik C., Erlangga Y.A. Comparison of Two-level Preconditioners Derived from Deflation, Domain Decomposition and Multigrid Methods J. Sci. Comput. 2009. V. 39. P. 340-370.
18. Domain Decomposition Methods.
URL: http://ddm.org
19. Il’in V.P. Metody i tehnologii konechnyh jelementov [Finite Element Methods and Technologies] Novosibirsk, ICMMG SB RAS, 2007.
20. Il’in V.P. Ob iteracionnom metode Kachmazha i ego obobwenijah [On Iterational Kachmazh Method and its Generalizations]. Sib. Zh. Ind. Mat. [J. Appl. Industr. Math.], 2006. V. 9, No 3. P. 39-49.
Поступила в редакцию 14 июня 2012 г.
