Метод Капорина-Коньшина параллельной реализации блочных предобусловливателей для несимметричных матриц в задачах фильтрации многокомпонентной смеси в пористой среде Текст научной статьи по специальности «Математика»

УДК 519.6

МЕТОД КАПОРИНА-КОНЬШИНА ПАРАЛЛЕЛЬНОЙ РЕАЛИЗАЦИИ БЛОЧНЫХ ПРЕДОБУСЛОВЛИВАТЕЛЕЙ ДЛЯ НЕСИММЕТРИЧНЫХ МАТРИЦ В ЗАДАЧАХ ФИЛЬТРАЦИИ МНОГОКОМПОНЕНТНОЙ СМЕСИ В ПОРИСТОЙ СРЕДЕ

Рассмотрены блочные предобусловливатели класса ILU (ILU(0), ILU(1), ILUT) для итерационных методов решения систем с разреженными матрицами, возникающими при аппроксимации систем дифференциальных уравнений в частных производных, описывающих фильтрацию многокомпонентной смеси в пористой среде. Предложена параллельная реализация блочного варианта построения ILU-разложения с хорошими характеристиками сходимости, имеющая значительное ускорение по сравнению с последовательной версией. Проведены численные эксперименты с использованием различных матриц, полученных при дискретизации реальных задач на моделях нефтяных месторождений Западной Сибири.

Ключевые слова: параллельное блочное ILU-разложение, разреженные матрицы, итерационный алгоритм, предобусловливатель.

The ILU class preconditioners (ILU(0), ILU(1), ILUT) employed for iterative algorithms for nonsymmetrical linear sparse matrix systems are considered. Test matrices used in this study are originated from discretization of systems of partial differential equations describing multicomponent fluid flow in porous media. Novel parallel algorithm for block ILU factorization is suggested. This algorithm demonstrates a good convergence and significant speed-up in comparison with sequential algorithms. New integrated approach was tested on the wide range of matrices resulted from real hydrodynamic simulations of oil fields of Western Siberia and demonstrated significant reduction in computational time.

Key words: parallel block ILU factorization, sparse matrices, iterative algorithm, preconditioner.

Введение. Для решения систем линейных алгебраических уравнений (СЛАУ) большой размерности на параллельных ЭВМ применяются различные методы разбиения матрицы между процессами и различные предобусловливатели. Классическим методом разбиения является блочно-диагональный [1], с помощью которого естественным образом распараллеливаются стандартные алгоритмы решения СЛАУ. И. Е. Капорин и И. Н. Коньшин [2] при решении симметричных положительно-определенных СЛАУ методом сопряженных градиентов предложили использовать для построения предобусловливателя метод перекрывающихся разбиений матрицы СЛАУ на блоки. К матрице применяется блочная версия двустороннего неполного обратного разложения Холецкого, а затем каждый из блоков заменяется на аппроксимацию — результат неполного разложения Холецкого. Полученный таким образом метод обладает лучшими характеристиками сходимости и параллелизма по сравнению с блочно-диагональным методом разбиения.

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

1. Представление предобусловливатели 1Ьи в аддитивной форме. Рассмотрим матрицу А

вида

К. Ю. Богачев1, Я. В. Жабицкий

2

1 Богачев Кирилл Юрьевич — канд. физ.-мат. наук, доцент каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: [email protected]. msu. su.

2Жабицкий Яков Вячеславович — асп. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: [email protected].

Пусть Ь и и — матрицы, соответствующие точному ЬИ-разложению матрицы А. Обозначим: Ь' = Ь 1, и' = и-1. Следуя [2], рассмотрим произведение

Н

и 1 и12 0 о и2 2 и2 з

оо из

зз

,ь\1 о о Ь21 Ь22 о о Ь32 Ь33/

1

А11 А12 -1 о

= А21 А22 о

о о

'о о о + I о (А22 А23

о

о V А32 Азз

-1

'о о о^ - I о А—1 о ко о о,

Введем обозначения:

А'' =

А22 А2з Аз2 Азз

а Ь'' и и'' — матрицы, соответствующие точному ЬИ-разложению матрицы А". Тогда

А22 А2з Аз2 Азз

1

2-21 о

о о

и'' и''

и22 и23

о и33

1

I о о I

Ь'2'2 о Ь'' Ь''

Ь32 Ь33

1

и'' и''

и22 и23

о и33

1

I о оо

Ь'2'2 о Ь'' Ь''

Ь32 Ь33

1

и'' и''

и22 и23 о и3''3

1

оо о I

Ь'2'2 о Ь'' Ь''

Ь32 Ь33

1

Таким образом, исходная матрица Н представлена в виде суммы двух матриц меньшей размерности, для каждой из которых можно независимо находить ЬИ-разложение.

Приведенные выкладки показывают, что матрица Н близка к матрице А-1 и является хорошим предобусловливателем для итерационных алгоритмов решения СЛАУ.

Заменив в выкладках точное ЬИ-разложение на 1ЬИ-разложение, мы получим аналогичный результат и для 1ЬИ-разложения.

Отметим также, что аналогично можно представить любую блочную трехдиагональную матрицу размерности (к = 4, 5,...) в виде суммы к — 1 матриц меньшей размерности, для каждой из которых можно независимо находить ЬИ (или 1ЬИ)-разложение.

2. Построение ХШ-разложения. Рассмотрим произвольную матрицу А размерности п. Пусть р — число процессов, П1,П2,... ,пр+1 — такие числа, что щ + П2 + ... + пр+1 = п. Очевидно, что, обнулив некоторые элементы матрицы А, можно получить блочную трехдиагональную матрицу А', такую, что ее г-й диагональный блок будет иметь размерность п^. Как было показано в п. 1, предобусловливатель для упомянутой блочной трехдиагональной матрицы А' можно представить в виде суммы р матриц меньшей размерности, для каждой из которых можно независимо найти 1ЬИ-разложение. Здесь мы не будем останавливаться на описании алгоритмов построения 1ЬИ-разложения, поскольку они подробно описаны в работах [1] (для скалярных версий) и [3] (для блочных версий).

3. Методы разбиения матрицы на части. Для дальнейшего изложения важно следующее замечание: чтобы вычислить элементы 1ЬИ-разложения, находящиеся в г-й строке, нам необходимо использовать уже рассчитанные элементы всех строк с меньшими номерами. При аддитивном представлении матрицы А в каждом слагаемом Ак (к = 1,...,р) есть строки, начинающиеся лишь с некоторого номера. Соответственно получаемые при 1ЬИ-разложении матриц Ак матрицы Ьк и ик при всех к > 1 значительно отличаются от соответствующих частей 1ЬИ-разложения исходной матрицы А.

Отметим, однако, что в силу разреженной структуры матриц при вычислении г-й строки 1ЬИ-разло-жения нужны только строки с номерами ненулевых элементов г-й строки.

Отсюда следует, что, с одной стороны, чем больше строк будет включено в каждую из подматриц, тем точнее получится 1ЬИ-разложение и соответственно будет меньшее число итераций, а с другой — чем больше строк мы добавим, тем существеннее увеличится объем вычислений и каждая итерация будет занимать большее время. Это означает, что существует оптимальное разбиение исходной матрицы на части.

3.1. Геометрический подход. Пусть задача решается в трехмерной области, аппроксимированной трехмерной сеткой, размерность которой по х, у и г составляет щ, п2 и п3 соответственно. Положим

I = а шш ТТ

¿€{1,2,3} -1-1

зе{1,2,з}

т.е. Ь — минимальное из попарных произведений т, увеличенное в а раз. Разобьем матрицу А на р равных частей и добавим Ь предыдущих строк к каждой части начиная со второй.

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

3.2. Алгебраический подход. У геометрического подхода есть существенный недостаток: в случае, когда физическая модель обладает сложной структурой, т.е. имеются связи между элементами, которые лежат не в соседних слоях сетки (например, сильно неравномерная сетка), данный подход может привести к получению плохого 1ЬИ-разложения, так как игнорируются существенные элементы исходной матрицы.

Алгебраический подход устраняет этот недостаток геометрической реализации.

Итак, пусть Ак — к-я часть матрицы А при равномерном разбиении, 10 — номера всех строк, попавших в Ак. Проведем первый шаг алгоритма. Обозначим . = 10 и € 10 : а^ = 0}, т.е. мы присоединили к уже имеющимся в множестве 1к индексам все номера столбцов ненулевых элементов, содержащихся в строках матрицы Ак. Матрицу, состоящую из строк с номерами из множества 11, обозначим А к. На втором шаге алгоритма описанную для множества 10 и матрицы Ак процедуру проведем для множества 1к и матрицы Ак, в результате получим множество 1к и матрицу Ак, и т.д. Таким образом можно построить последовательность множеств Ц}?=1 и матриц [Агк}?=1. Очевидно, что в силу конечности строк исходной матрицы А начиная с некоторого I все члены указанных последовательностей совпадут.

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

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

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

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

Рис. 2. Пример распределения между процессами элементов сетки для прямоугольной области с вырезом: а — равномерное разбиение; б — геометрический подход при а = 1; в — алгебраический подход (1 шаг) (серым цветом отмечены строки, присоединенные ко второму процессу, а клеткой — к третьему)

4. Численные эксперименты. Численные эксперименты с применением описанных выше алгоритмов ILU, ILU(1) и ILUT(p, т) проводились на матрицах, получаемых при дискретизации реальных задач течения жидкости (использовались углеводородные месторождения Западной Сибири) на сетках большой размерности как регулярной, так и нерегулярной структуры. Задачи запускались на системе с двумя процессорами Intel Nehalem (2xXeon 5570) c общей памятью 20 Гб, частотой каждого процессора 3,0 ГГц и на системе "Regatta" Московского государственного университета (IBM eServer pSeries 690) с общей памятью 64 Гб, состоящей из 16 процессоров Power4, частота каждого процессора 1,3 ГГц. На системе Intel Nehalem задачи запускались на 1, 2, 4 и 8 процессорах, а на системе "Regatta" — на 1, 2, 4, 8 и 15 процессорах (один процессор занимает операционная система, поэтому запуск на 16 ядрах не дает ускорения, а лишь искажает результаты численных экспериментов).

Постановка задачи. Рассматриваемые задачи сводятся к следующей системе уравнений [4]:

О т

- (фКс) = div .re,/£z(k^(V№-7zVI))) +qc,

Po Pg - Pc og , Po Pw - Pc ow j C - 1, ... , nc ,

Sw + So + Sg — 1,

Bw V Bo BgJ \Bg Bo

(1)

где

nc — количество компонентов смеси;

Rg,o — Rg,o(po) — растворимость газа в нефтяной фазе;

Bl — Bi(pi) — коэффициент объемного расширения фазы l, l — o (oil), g (gas), w (water); Ф — Ф(PwJ Po, Pg, x, y, z) — пористость среды;

xc,l — xCll(pWJp0JPgJ N) — молярная доля компонента c в фазе l; N — (Ni,Nnc) — вектор молярных плотностей; Cl — Cl(Pl, N) — молярная плотность фазы l; k — k(pWJp0JPgJЖJ y,z) — тензор абсолютной проницаемости; krl — krl(Sw,Sg) — относительная проницаемость фазы l; ßl — ßl(pi) — вязкость фазы l;

Yl — Pl9 — вертикальный градиент гидростатического давления в фазе l; g — гравитационная постоянная; D — D(x,y,z) — вектор глубины (сверху вниз); Pl — Pl(Pl) — массовая плотность фазы l;

Pcog — Pcog(Sg) — капиллярное давление в системе нефть-газ; Pcow — Pcow(Sw) — капиллярное давление в системе вода-нефть; qc — qc(pw,po,Pg, NJÍJXJУJz) — источник компонента c (скважина). Неизвестными в этой системе являются:

1) Nc — Щ(Ь,х,у, z) — молярная плотность компонента c (для модели черной нефти компонентами служат вода, нефть и газ);

2) Sl — Sl(t, х, у, z) — насыщенность фазы l;

3) Pl — Pl(t, х, у, z) — давление в фазе l.

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

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

F(p, Ni,...,Nnc) —0,

где p — (рг), Nc — (NO — векторы значений давления и молярных плотностей в блоках сетки. Затем для решения системы нелинейных уравнений F(U) — 0, U = (p, N) применяется стандартный метод Ньютона:

ljm+i = иш_ ^Щр^ F(Um).

Здесь dF(Um)/dU — отображение (матрица) Rnc*(^+J) ^ Rnc*(^+J) х Rnc*(^+J) (одна переменная исключается с помощью третьего уравнения (1)), K — число блоков сетки, J — число скважин. На каждом шаге метода Ньютона необходимо решать систему с несимметричной матрицей dF(Um)/dU, т.е. задача сводится к решению системы линейных уравнений

Ax — r, (2)

матрица которой есть якобиан метода Ньютона. Матрицу A можно рассматривать как матрицу, элементами которой являются блоки размера (nc) х (nc). Значение nc в свою очередь варьируется в зависимости от типа задачи от 2 до 21.

Для решения системы линейных уравнений (2) используется итерационный алгоритм BCGSTAB [1], в котором в качестве предобусловливателя берутся блочные варианты [3] рассматриваемых ILU-разложений: ILU-разложение с использованием геометрического подхода при а — 2 (ILUG) и алгебраического подхода с двумя шагами (ILUA); ILU^-разложение; ILUT-разложение (с параметрами p — 5 и

Т =

10 3). В качестве критерия окончания итерационного процесса выбрано условие падения относитель-

ной невязки на 6 порядков.

Для численных экспериментов взяты матрицы системы линейных уравнений (2) дискретизации реальных физических моделей (1) на разных шагах по времени.

полученные при

Характеристики матриц

№ пх Пу nz d г Ъ п то ksk

1 120 160 6 2 199490 649075 2596300 И 3,564 • 105

2 83 107 85 2 205650 699501 2798004 49 2,234 • 107

3 126 109 54 2 283996 801368 3205472 12 1,146 • 108

4 130 123 17 2 382580 1200506 4802024 18 1,226 • 106

5 114 238 12 2 364002 1200577 4802308 15 7,955 • 107

6 88 71 28 3 281457 552287 4970583 24 1,620 • 107

7 133 120 50 2 499322 1491985 5967940 47 1,813 • 107

8 133 120 50 2 568006 1669725 6678900 36 4,108 • 105

9 392 495 3 2 565086 1786979 7147916 7 1,340 • 105

10 71 71 575 2 759746 1905443 7621772 7 1,632 • 106

И 194 121 55 2 589580 1905520 7622080 40 4,786 • 106

12 122 147 35 2 679262 2190961 8763844 22 4,276 • 106

13 138 81 50 2 748776 2437818 9751272 40 1,395 • 108

14 46 99 173 3 510750 1156672 10410048 14 1,860 • 1011

15 97 135 30 3 575025 1205459 10849131 7 2,139 • 103

16 133 103 117 2 815932 2743130 10972520 53 3,042 • 105

17 80 71 125 2 858382 2943033 11772132 63 2,904 • 106

18 398 204 3 3 730728 1539038 13851342 7 2,550 • 103

19 90 174 42 3 702405 1609689 14487201 7 9,891 • 105

20 60 55 276 3 766098 1721412 15492708 17 1,821 • 106

21 276 381 12 2 1391786 4242249 16968996 14 2,362 • 106

22 420 280 30 2 1353710 4536033 18144132 30 1,396 • 106

23 99 165 126 3 895212 2043612 18392508 10 4,724 • 105

24 181 116 83 3 1028625 2240069 20160621 39 1,024 • 107

25 78 263 31 3 1121196 2456774 22110966 23 7,889 • 109

26 389 244 30 2 1749732 5726484 22905936 28 1,400 • 108

27 249 198 12 3 1393869 2980557 26825013 16 3,651 • 106

28 190 134 575 2 2702850 6746339 26985356 7 2,847 • 103

29 407 412 37 2 2400402 7367991 29471964 28 1,268 • 106

30 60 220 85 2 2188842 7484949 29939796 88 2,486 • 106

31 219 454 150 2 2715768 8579656 34318624 65 3,658 • 107

32 117 214 64 3 1965948 4055996 36503964 45 3,096 • 101и

33 354 750 36 3 2348076 4921354 44292186 10 6,782 • 101и

34 198 366 104 2 4863578 15818315 63273260 48 5,120 • 106

35 280 290 80 3 4970010 8818554 79366986 18 5,144 • 106

36 88 215 177 3 7256967 16524695 148722255 9 5,475 • 106

Результаты. В таблице приведены некоторые характеристики матриц: пх — размер сетки по х; пу — размер сетки по у; пх — размер сетки по г; I — размер блоков; г — количество скалярных строк матрицы (соответственно количество блочных строк равно г/1); Ь — количество хранящихся блоков; п — количество хранимых скалярных элементов (равно Ь ■ I2); т — максимальное количество ненулевых блоков в строке матрицы; к3и — коэффициент кососимметричности, вычисленный по формуле

ksk = max

i=l,...,r/d

r/d

i]

— A

=i

где || ■ — строчная норма матрицы.

Пусть — время работы (или количество итераций) j-го алгоритма на р процессах на 1-й матрице. Тогда рассмотрим следующее усреднение:

N

JiPi = J2P2

1 Ñ

£

i=1

а

Jipi

J2P2

а

(3)

здесь Pi,P2 — число процессов, на которых запускался расчет задачи, а ji ,j2 обозначают алгоритм, использующийся при решении, и принимают значения ILUG, ILUA, ILU(1), ILUT. Данный коэффициент cj^^l показывает среднее ускорение (или среднее увеличение числа итераций).

На рис. 3, 4 приведены результаты сравнения между собой алгоритмов с использованием ILUG-, ILUA-, ILU(1)- и ILUT(5, 10_3)-разложений при расчете на системе Intel Nehalem, а на рис. 5,6 — сравнение тех же алгоритмов при расчете на системе "Regatta".

Как показывают проведенные численные эксперименты, средний рост числа итераций не превышает 6% для всех рассматриваемых алгоритмов, что говорит о том, что благодаря использованию метода Капорина-Коньшина параллельной реализации блочных предобусловливатей получаемое ILU-разложение обладает хорошими характеристиками масштабируемости. Наибольший рост числа итераций не превышет 40% при расчете на 8 ядрах.

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

Отметим также, что за счет расширения шаблона ненулевых элементов число итераций для ILU(1)-и ILUT-разложений существенно меньше, чем для ILUG- и ILUA-разложений.

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

Рис. 3. Сравнение среднего (усреднение бралось по формуле (3)) количества итераций при расчете на системе Intel Nehalem со следующими значениями параметров ji, pi. j = ji € {ILUG, ILUA, ILU(1), ILUT}, pi € {1;2;4;8}, P2 = 1 (а); ji = ILUG, j2 € {ILUG, ILUA, ILU(1), ILUT}, pi € {1; 2; 4; 8}, p2 = 1 (б)

Рис. 5. Сравнение среднего (усреднение бралось по формуле (3)) количества итераций при расчете на системе "Regatta" со следующими значениями параметров ji, pi. j2 = ji € {ILUG, ILUA, ILU(1), ILUT}, pi €{1; 2; 4; 8; 15} p2 = 1 (а); ji = ILUG, j2 € {ILUG, ILUA, ILU(1), ILUT}, pi € {1; 2; 4; 8; 15}, p2 = 1 (б)

Рис. 4. Сравнение среднего (усреднение бралось по формуле (3)) ускорения при расчете на системе Intel Nehalem со следующими значениями параметров ji, pi. ji = j2 € {ILUG, ILUA, ILU(1), ILUT}, pi = 1, p2 € {1; 2; 4; 8} (а); ji € {ILUG, ILUA, ILU(1), ILUT}, j2 = ILUG, pi = 1, p2 € {1; 2; 4; 8} (б)

Рис. 6. Сравнение среднего (усреднение бралось по формуле (3)) ускорения при расчете на системе "Regatta" со следующими значениями параметров ji, pi. ji = j2 € {ILUG, ILUA, ILU(1), ILUT}, pi = 1, p2 €{1; 2; 4; 8; 15} (а); ji € {ILUG, ILUA, ILU(1), ILUT}, j2 = ILUG, pi = 1, p2 € {1; 2; 4; 8; 15} (б)

метим здесь, что при использовании 1ЬИ(1)-разложения характиристики параллелизации хуже, чем при использовании остальных алгоритмов: это объясняется тем, что стратегия расширения шаблона ненулевых элементов предобусловливателя никак не учитывает величины добавляемых элементов, что приводит к увеличению арифметических операций, но не к повышению точности разложения. ILUT-разложение устраняет данный недостаток ILU-разложения и как следствие обладает лучшими характеристиками па-раллелизации.

Отметим, что на одном процессе в среднем ILUG- и ILUA-разложения гораздо предпочтительнее, чем ILU(1)- и ILUT-разложения. В первую очередь это объясняется вычислительными затратами на построение предобусловливателя и большим числом ненулевых элементов в предобусловливателе для ILU(1)- и ILUT-разложений, в результате чего каждая итерация занимает большее время. Однако с увеличением количества процессов картина заметно меняется: качество предобусловливателя играет все более значимую роль и время, затрачиваемое на итерационный процесс, компенсирует время, затраченное на построение предобусловливателя. При увеличении числа расчетных ядер в среднем наиболее предпочтительным оказывается ILUA-разложение.

СПИСОК ЛИТЕРАТУРЫ

1. Saad Y. Iterative Methods for Sparse Linear Systems. Philadelphia: SIAM, 2003.

2. Капорин И.Е., Коньшин И.Н. Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки // Журн. вычисл. матем. и матем. физ. 2001. 41, № 4. 515-528.

3. Богачев К.Ю., Жабицкий Я.В. Блочные предобусловливатели класса ILU для задач фильтрации многокомпонентной смеси в пористой среде // Вестн. Моск. ун-та. Матем. Механ. 2009. № 5. 19-25.

4. Aziz K., Settari A. Petroleum Reservoir Simulation. London: Applied Science Publishers, 1979.

Поступила в редакцию 27.05.2009