Алгебраический многоуровневый метод AMG: сравнение с методом BiCGStab + ILU и использование в составе метода CPR Текст научной статьи по специальности «Математика»

УДК 519.6

АЛГЕБРАИЧЕСКИЙ МНОГОУРОВНЕВЫЙ МЕТОД AMG: СРАВНЕНИЕ С МЕТОДОМ BICGSTAB + ILU И ИСПОЛЬЗОВАНИЕ В СОСТАВЕ МЕТОДА CPR

К.Ю. Богачев1, М. Ю. Михалева2, И. Г. Горелов3

Проведено сравнение на модельной задаче с сильной анизотропией алгебраического многоуровневого метода AMG и метода BICGSTAB+ILU. Метод AMG был использован в составе метода CPR для решения задачи фильтрации вязкой сжимаемой жидкости в пористой среде. Алгоритм CPR+AMG на данной задаче сравнивался с методом CPR+ILU. Проводился анализ работы CPR+AMG при решении задачи для модельного месторождения.

Ключевые слова: алгебраический многоуровневый метод, метод AMG, метод CPR, SPE9, течение вязкой жидкости в пористой среде.

Algebraic multigrid method AMG is compared with BICGSTAB+ILU method concerning to the model problems with strong anisotropy. AMG is considered in conjunction with CPR method concerning to the filtration problems of multicomponent fluid flow in porous media. CPR+AMG method is compared with CPR+ILU. CPR+AMG is considered in detail concerning to the model oil field.

Key words: algebraic miltigrid, AMG method, CPR method, SPE9, fluid flow in porous media.

Введение. Алгебраический многоуровневый метод AMG (Algebraic Multigrid) является одним из наиболее популярных методов решения систем линейных уравнений, полученных при аппроксимации дифференциальных уравнений методами конечных элементов или конечных разностей [1, 2]. Этот метод и его модификации были детально рассмотрены в статьях [3, 4]. Авторами этих статей был реализован алгоритм на языке FORTRAN (код AMG1R1).

В настоящей работе описан метод AMG и указываются места его модификации для примения в составе метода CPR. Алгоритм AMG реализован на языке C+ + , и проведено сравнение с методом би-сопряженных градиентов с предобусловливателем — неполным LU — на примере задачи, моделирующей распределение тепла в сильноанизотропной среде. Кроме того, исследовано качественное поведение алгоритма в зависимости от степени анизотропии.

Метод AMG используется в составе метода CPR (Contrained Pressure Residual [5]) для решения задачи фильтрации вязкой сжимаемой жидкости в пористой среде [6]. Комбинация методов CPR + AMG выступает в качестве предобусловливателя для системы линейных уравнений, получаемой на каждом шаге метода Ньютона решения задачи фильтрации по неявной схеме. Метод CPR + AMG на примере данной задачи сравнивается с методом CPR + ILU. Проводится анализ работы CPR + AMG при решении задачи для модельного месторождения при различных значениях параметров модели.

Описание алгоритма AMG. Подробное описание алгоритма приведено в статьях [3, 4]. Рассматривается система линейных уравнений A * u = f. Для применения алгебраического многоуровневого метода матрица A должна удовлетворять следующим условиям:

A — положительно-определенная квадратная матрица;

A имеет "существенно" положительный тип (essentially positive type),

т.е.

A — разреженная матрица;

все диагональные элементы матрицы положительны;

большая часть внедиагональных элементов меньше нуля или равна нулю;

суммы элементов в строке больше нуля или равны нулю (слабое диагональное преобладание).

1 Богачев Кирилл Юрьевич — канд. физ.-мат. наук, доцент каф. вычислительной математики мех.-мат. ф-та МГУ^-mail: bogachev@mech.math. msu. su.

2Михалева Майя Юрьевна — асп. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: maya-mm@yandex.ru.

3Горелов Илья Геннадьевич — асп. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: chronosphere@mail.ru.

Рассмотрим двухуровневый алгебраический алгоритм [3, 4] решения системы линейных уравнений A * u = f (аналогично строится многоуровневый).

В алгебраическом методе используется терминология фиктивных сеток (h и H — мелкая и грубая сетки). Пусть i-я точка — это каждый элемент i Е = {1, 2,... ,n}, i — номер строки (столбца)

исходной матрицы. Матрица A представляется в виде графа связей точек. Точка i Е связана с точкой j Е Qh, если ahj = 0. Систему линейных уравнений можно записать как сеточное уравнение на фиктивной

мелкой сетке Qh: A^uh = fh.

На грубой сетке AhuH = fH, где Ah = IHAhIj.

Двухуровневый алгоритм имеет вид uhew = u^ + Ij eH, где Ah eH = rH = Ij (r£ld) = Ij (fh — AhU^).

Метод AMG можно разделить на две фазы — инициализация и решение:

фаза инициализации включает в себя построение F/C-разбиения и операторов Ij и Ij для каждой сетки. При этом разбиение и операторы должны быть такими, чтобы оператор Ah оставался разреженным и был значительно меньшей размерности, чем Ah, и чтобы достигалась достаточная сходимость алгоритма. Фаза инициализации может быть реализована один раз — в случае, когда AMG используется в качестве предобусловливателя в итерационном алгоритме;

фаза решения включает в себя итерационный процесс, каждая итерация которого представляет собой V-, W- или F-цикл, сглаживание методом Гаусса-Зейделя и решение системы уравнений на самом грубом уровне.

Сравнение методов AMG и BICGSTAB+ILU на задаче с сильной анизотропией. Описанный алгоритм был реализован на языке C++, и проведены численные эксперименты на задаче с сильной анизотропией.

В табл. 1 и 2 представлены результаты работы программы в сравнении с методом BICGSTAB+ILU [7]. Испытания проводились на компьютере AMD Atlon 1833 МГц, 512 кбт.

Рассмотрим задачу —div(kVu) + \u = f с нулевыми краевыми условиями в области [0,1] х [0,1]. Будем использовать равномерную сетку, число узлов по x и по y равно 4096.

Пример 1. Возьмем матрицу k следующего вида: k = 1 1 0

0 w

Таблица!

Метод Параметр w

3 100 500 2000 10000 40000

AMG 2(10) 2(8) 1(7) 1(7) 1(7) 1(7)

BICGSTAB+ILU 8(53) 13(87) 23(148) 36(228) 61(373) 80(471)

При различных значениях параметра и> сравниваются скорости работы алгоритмов. В AMG используется V-цикл. В табл. 1 представлено сравнение времени работы и числа итераций методов (число итераций указано в скобках). Из таблицы видно, что с возрастанием параметра и> алгоритм AMG все больше и больше превосходит по скорости работы метод BICGSTAB + ^и.

Пример 2. Теперь рассмотрим случай W-цикла, к — матрица следующего вида: к =

Таблица2

Метод Параметр е

0,001 0,01 0,03 0,035 0,038 0,042

AMG 4(10) 5(10) 7(16) 11(24) 20(40) Расходится

BICGSTAB+ILU 10(69) 10(68) 12(71) 12(81) 13(85) 21(138)

При малых е метод AMG превосходит по скорости BICGSTAB + ^и, однако при увеличении параметра е он постепенно теряет свои преимущества и начинает проигрывать (табл. 2). Далее алгоритм AMG начинает расходиться, тогда как BICGSTAB + ^и продолжает сходиться.

Применение ЛИС в составе метода CPR для решения задачи фильтрации. Задача фильтрации вязкой сжимаемой смеси в пористой среде (подробно постановка задачи фильтрации приведена в работе [6]) записывается следующим образом:

д к

— (фМс)=(Иу V" хс>р^р13(\и-^{У'Рр - ррдУО)) + qc, с = 1,...,пс, пс = 3, дг \ ар ) /1\

р=о,ш,о гр (1)

РО - Ро = РеОС, Ро - Рш = РеОШ,

13 ВМУ, математика, механика, №4

Бп? + ¿>о + Бо = 1, N1 = =

М2 = Мо= &М1Г + М3 = мс= &М1Г +

Во Во Во Во

Здесь искомыми являются следующие величины:

1) Нс = МС(Ь, х, у, г) — молярные плотности;

2) Рр = Ррх, у, г) — давление (Р = 0,Ш,0 — нефтяная, водная и газовая фазы);

3) Бр = Бр(Ь,х,у,г) — насыщенности фаз.

Величины ф, хс,р, , к, кГрр, ¡р, рр, дс, РСоо, РСош, Вр, Ко,о, Но,о являются заданными функциями, нелинейно зависящими от искомых величин, и определены в (1).

Обычно первое уравнение (с = 1) заменяется на сумму уравнений по всем с, а именно на уравнение (при суммировании используется то, что ^ хс,р = 1)

с

| (фЕЛс) = <Иу £ ~ РРдУО)) + (2)

01 \ с ) Р=о,Щ,С ¡Р С

С помощью неявной схемы (для аппроксимации по времени) и метода конечных объемов (для аппроксимации по пространственным переменным) задача сводится к системе нелинейных уравнений

Е(р, ) = 0,

где р = (рг), N = (N0 — векторы значений давления и молярных плотностей в блоках сетки. Для решения системы используется метод Ньютона

ит+1 = иш_ ^Шр^ Р{ит),

где Г(и) =0, и = (р, N); дГ(ит)/ди — матрица К"с*(к+1) ^ кпс*(к+^ х Кпс*(к+^, к — число блоков сетки, .] — число скважин. На каждом шаге метода Ньютона решается система с матрицей дГ(ит)/ди, задача сводится к решению системы линейных уравнений Ах = г, А — якобиан из метода Ньютона. Матрицу А можно рассматривать как матрицу, элементами которой являются блоки размера (пс) х (пс). Значение пс в свою очередь варьируется в зависимости от типа задачи от 2 до 21.

Далее применяется предобусловливание методом СРИ (глобальный алгоритм) [5]. В силу того что наибольший вклад в вектор невязки вносит часть матрицы, соответствующая переменной давления, во всех внедиагональных блоках приравниваются к нулю производные по переменным, отличным от р.

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

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

Обозначим процесс коррекции методом Гаусса через Gauss(A), а процесс выделения из матрицы А подматрицы давления — через Иеёисе(А). Тогда матрица для приближенного вычисления давления Ар находится по формуле Ар = Кеёисе^аи88^и11(А))).

Для решения подсистемы с матрицей Ар применяется метод AMG (т.е. используется комбинация методов СРИ + AMG). Практическая реализация такого подхода с некоторыми упрощениями рассмотрена в работе [8].

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

Применимость AMG зависит от "почти эллиптичности" системы по переменной давления. По закону Дарси поток фазы

кгр

иР = ~(Зс(к—{УрР - ррдУО)).

V ¡¡Р }

Уравнение (2), записанное через поток фазы, примет вид

• + div £ №> (3)

с \ с ) с P=O,W,G

где С = — — сжимаемость породы. Согласно данной записи, при отсутствии скважин (параметр qc), dt

ускорения свободного падения (g) и капиллярного давления (Pcow,Pcog) уравнение (3) запишется следующим образом:

cl+4 (?"")=div(aVp)l

a = a(),,N)= (,3 £ ip^V,

У P=O,W,G VP J

где a — симметричная, положительно-определенная матрица, что обеспечивает эллиптичность задачи для давления при неявной аппроксимации по p и явной по Nc.

Зависимость задачи от каждого из параметров задачи показана на примере следующего численного эксперимента. Модель — синтетический трехфазный тест небольшого размера SPE9, предлагаемый Сообществом инженеров-нефтяников (Society of Petroleum Engineers, или SPE) как эталонный для тестирования работоспособности программ гидродинамического моделирования.

Расчет производится для 899 дней. "Отчетные" временные шаги достаточно большие, при расчете разделяются на более мелкие — "расчетные". Длины "отчетных" шагов представлены в табл. 3.

ТаблицаЗ

Длина, день Номер шага

1 2 3 4 5 6 7 8 9 10 И 12 13 14 15

1 1 1 47 50 100 100 60 60 60 60 60 60 60 179

Расчет стартует с состояния равновесия. Параметры будем варьировать следующим образом:

1) длину шага по времени ограничим значением 10 дней (ввиду того что величина расчетного шага выбирается программой автоматически, в частности не длиннее текущего "отчетного" шага, данное значение может и не достигаться);

2) ускорение свободного падения положим равным нулю;

3) скважины удалим;

4) капиллярное давление положим равным нулю.

Данные упрощения применяем в различных сочетаниях к модели и решаем ее методом СРИ. В качестве алгоритма решения локальной системы используем AMG или BCGStab+ILU(0). Результаты тестирования приведены в табл. 4.

Таблица4

Применяемые упрощения и количество шагов метода Номер варианта

0 1 2 3 4 5 8 9 10 И 12 13

Ограничение длины шага + + + + + +

Отключение гравитации - - - - + + - - - - + +

Удаление скважин - - + + - - - - + + - -

Выключение капиллярного давления - + - + - + - + - + - +

СРК+АМС, шаг 4 4 4 15 4 4 10 10 4 15 6 4

СР11+1Ш, шаг 15 15 15 15 15 15 15 15 15 15 15 15

Знак "+" указывает на наличие упрощения, а знак " —" — на отсутствие. В двух последних строкак приведено число пройденных отчетных шагов до первого случая расходимости. Если указано 15 отчетных шагов, расчет дошел до конца модели. Варианты 6, 7, 14, 15, соответствующие удалению скважин и отключению гравитации, являются тривиальными, поскольку для них решение задачи — это нулевое приближение в методе Ньютона.

14 ВМУ, математика, механика, №4

Расходимость во всех вариантах, где она имела место, вызвана отсутствием сходимости локального алгоритма, поскольку замена AMG на BCGSTAB+ILU(0) приводит к тому, что расчет доходит до конца модели (последняя строка таблицы).

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

Кроме того, когда упрощением было только ограничение по длине шага (вариант 8), алгоритм посчитал больше шагов (10), чем когда к нему добавлялось удаление скважин (вариант 10) — 4 шага или "отключение гравитации" (вариант 12) — 6 шагов. Это обусловлено тем, что усложнение матрицы приводит к тому, что метод автоматического вычисления шага не позволяет расчету сразу "разогнаться" до длины в 10 дней.

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

Заключение. 1. Метод AMG был протестирован на модельных задачах, близких к симметричным эллиптическим, и показал хорошие результаты расчета, причем тем лучше, чем более симметрична сама задача.

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

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

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

2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Лаборатория Базовых Знаний, 2003.

3. Stuben K. A Review of Algebraic Multigrid. GMD Report 69, November, 1999.

4. Stuben K. Algebraic Multigrid (AMG): an Introduction with Applications. GMD Report 70, 1999.

5. Wallis J.R. Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration. SPE 12265, 1983.

6. Азиз Х, Сэттари Э. Математическое моделирование пластовых систем / Пер. с англ. М.: Недра, 1982.

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

8. Lacroix S., Vassilevski Yu., Wheeler J., Wheeler M. Iterative solution methods for modeling multiphase flow in porous media fully implicitly // SIAM J. Sci. Comput. 2003. 25, N 3. 905-926.

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

УДК 519.21

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

Н. А. Ирхина1

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

Ключевые слова: принцип подсчета премии, функция искажения, абсолютная чувствительность премии.

This research is devoted to Wang's premium principle in actuarial theory. By example of Pareto distribution the author notes that Wang's premium principle can be applied to

1 Ирхина Наталья Александровна — асп. каф. теории вероятностей мех.-мат. ф-та МГУ, e-mail: irkhina.natalia@rambler.ru.