Численное исследование предобуславливателя Generalized Nested Factorization для задач пластовой фильтрации Текст научной статьи по специальности «Математика»

УДК 519.6

Численное исследование предобуславливателя Generalized Nested Factorization для задач пластовой фильтрации

В. Е. Борисов*, Е. Б. Савенков+

* Механико-математический факультет МГУ им. М.В. Ломоносова Ленинские горы, д. 1, Москва, 119991, Россия ^ Институт прикладной математики им. М.В. Келдыша РАН Миусская пл., д. 4а, Москва, 125047, Россия

Работа посвящена численному исследованию метода предобуславливания Generalized Nested Factorization (GNF) для задач пластовой фильтрации, отличительной особенностью которых является использование неструктурированных сеток и наличие нелокальных связей между ячейками.

Приведены результаты сравнительного анализа алгоритма со стандартными предобу-славливателями на основе неполного LU-разложения (ILU(0), ILU(1), ILUT). Численно исследованы спектральные свойства предобусловленных матриц.

Результатом работы является вывод о хорошей робастности алгоритма GNF для реальных задач пластовой фильтрации.

Ключевые слова: пластовая фильтрация, предобуславливание, неполное LU-раз-ложение, итерационные методы решения СЛАУ, generalized nested factorization.

1. Введение

В ходе разработки нефтегазового месторождения неизбежно возникает необходимость численного моделирования происходящих в нем процессов. Целью моделирования, например, может являться оптимизация процесса добычи или расчёт прогнозных вариантов разработки месторождения.

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

Основу рассматриваемого класса математических моделей составляют нелинейные уравнения фильтрации многофазного многокомпонентного флюида в пористой среде. Использование метода Ньютона для решения соответствующей нелинейной системы разностных уравнений приводит к необходимости многократного решения систем линейных алгебраических уравнений (СЛАУ), возникающих на каждой итерации метода Ньютона на каждом временном шаге. В силу ряда особенностей задач пластовой фильтрации (аспектные отношения ячеек и скачки коэффициентов ~ 10-100, наличие в системе уравнений, которые описывают сильно различающиеся по скорости процессы и др.), возникающие СЛАУ являются плохо обусловленными системами общего вида (несимметричными и не положительно определёнными) и обладают большой размерностью (до ~ 106 — 109 уравнений). Эти причины затрудняют их решение, существенно увеличивая время расчёта. Поэтому эффективность решение СЛАУ является критичной для эффективной работы программы-симулятора.

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

Статья поступила в редакцию 30 октября 2012 г.

Работа выполнена при частичной финансовой поддержке Минобрнауки России (госконтракт № 07.514.11.4008) и РФФИ (проект № 12-01-00793-а).

Авторы хотели бы выразить свою признательность к.ф.-м.н. Д.Ю. Максимову (ИПМ РАН) за советы и полезные обсуждения.

Эта работа не была бы выполнена без поддержки и внимания к.ф.-м.н. А.Х. Пергамент (ИПМ РАН).

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

Рассмотрим СЛАУ Ах = b с невырожденной квадратной матрицей А £ Rnxn и правой частью b £ Rn. Пусть В £ — невырожденная матрица. Тогда

система В-1 Ах = В-1Ь будет алгебраически эквивалентна исходной при любом выборе матрицы В.

Смысл предобуславливания заключается в выборе такой матрицы В, что

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

2. Решение системы уравнений By = z вычисляется существенно быстрее, чем решение исходной системы.

Матрица В называется матрицей предобуславливания, или просто предобуслав-ливателем [2].

В настоящее время известно много разных типов предобуславливателей, например, предобуславливатели на основе аппроксимации матрицы системы: семейства ILU, IQR и ILQ; предобуславливатели на основе аппроксимации обратной матрицы: полиномиальные, разреженные аппроксимации обратной матрицы (например, AINV), аппроксимации обратной матрицы в факторизованной форме (такие как FSAI, SPAI и др.) [3].

При решении задач пластовой фильтрации наиболее часто используются пре-добуславливатели семейства ILU, в частности, различные модификации алгоритма Nested Factorization, а также иерархические предобуславливатели типа CPR (Constrained Pressure Residual) [4].

Предобуславливатель Nested Factorization (NF) — эффективный предобуслав-ливатель, являющийся разновидностью блочного ILU-предобуславливания [5]. Его отличительной особенностью является малое время построения и умеренное количество используемой памяти. При построении предобуславливателя NF используется тот факт, что при аппроксимации задач на регулярных сетках с помощью стандартных двухточечных аппроксимаций потока через грани ячеек, матрица системы уравнений имеет рекурсивно блочный трёхдиагональный вид. При использовании структурированных сеток и стандартных аппроксимаций потоков предобуславливатель NF очень эффективен.

Цель настоящей работы — программная реализация и численное исследование предобуславливателя Generalized Nested Factorization (GNF) [6], применимого в случае использования неструктурированных сеток и при наличии нелокальных связей между ячейками сетки. В настоящей работе продемонстрирована его численная эффективность при расчёте ряда тестовых и реальных задач пластовой фильтрации. Отдельное внимание уделено вопросу влияния рассмотренных предобуславливателей на спектральные свойства предобусловленных матриц. Несмотря на то, что предобуславливатель GNF был предложен сравнительно давно, результаты численного исследования его эффективности и применимости при решении практических задач моделирования пластовых систем авторам не известны.

2. Предобуславливатель Nested Factorization

Предобуславливатель Nested Factorization (NF) предложен в работе [5]. Он представляет собой разновидность блочного ILU-предобуславливания, с тем отличием, что матрица предобуславливателя имеет рекурсивно трёхдиагональный

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

Рассмотрим прямоугольную сетку, состоящую из Nx х Ny х Nz расчётных ячеек. Зададим порядок обхода ячеек сетки, зафиксировав какую-либо перестановку (а, ) координат х, у и z: первой меняющейся координатой станет координата вдоль оси а, затем последовательно вдоль осей ¡3 и 7 .В общем случае эффективность предобуславливания методом NF сильно зависит от выбора ориентации координатных осей. В качестве главной оси (по которой происходит обход ячеек сетки в последнюю очередь) выбирается та координатная ось, вдоль которой коэффициенты уравнений системы связаны наиболее сильно — в большинстве практических задач фильтрации это будет ось z. В соответствии с алгоритмом метода ячейки сетки объединятся в «столбцы» (0, 7 = const), которые, в свою очередь, объединяются в слои ячеек (7 = const). Совокупность всех слоев составляет полную сеточную область.

В этом случае матрицу А можно представить в виде

А = d + и + I + v + т + w + п,

где d — главная диагональ матрицы; и и I — верхняя и нижняя (блочные) диагонали, соответствующие связям ячеек в пределах одной «линии» ячеек, v и т — верхняя и нижняя (блочные) диагонали, соответствующие связям между столбцами ячеек в пределах одного слоя, и, наконец, w и п — верхняя и нижняя диагонали, соответствующие связям между слоями сетки.

Для сетки 4 х 4 х 2 структура матрицы А изображена на рис. 1.

Рис. 1. Сетка (слева) и структура матрицы В (справа) для сетки 4 х 4 х 2 при использовании предобуславливателя КЕ ( диагонали: «о» — п, «+» —

га, «х» — I, «*» — «0» — и, «0» — V, «0» — )

Предобуславливатель В ^ А методом NF строится в следующем виде:

В = (Р + п)(1 + Р-1 ь)), (1)

Р = (Т + т)(1 + Т-1 V), Т =(д + 1)(1 + д-1и), д = d — I д-1и — £, е = В — А = тТ-1ь + пР-1

(2)

(3)

(4)

ю,

где I — единичная матрица; Р — блочная трёхдиагональная матрица, блоки которой соответствуют «линиям» ячеек в пределах одного слоя; Т — трёхдиагональная матрица, элементы которой соответствуют связям отдельных ячеек в пределах отдельных «линий» ячеек; д — модифицированная главная диагональ матрицы А; е — матрица ошибки.

Отметим, что матрица В имеет блочную трёхдиагональную структуру, соответствующую структуре матрицы А.

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

В итоге предобуславливатель МР определяется соотношениями (1)—(4) при замене матрицы ошибки на её аппроксимацию . В этом случае получаем

В = А + £ — £.

Наиболее распространёнными являются следующие варианты выбора ё:

р(1) = р • р(2) = ь гг ь гг; ь гг

N

со18иш(е ) = ^^ е^;

,(3)

N

£■/ = ГО'№8ит(£) =

( )=

г]'

1=1

3 = 1

Первый способ даёт сравнительно грубое приближение и редко используется на практике.

Рассмотрим второй способ. Пусть X — решение системы Вх = Ь. Тогда невязка исходной системы на векторе х имеет вид: г = (В — А)Х. При этом, в силу построения е(2), всегда имеем со1виш(В — А) = 0, откуда следует, что сумма всех компонент вектора невязки равна нулю вне зависимости от значений х. Это свойство обеспечивает сохранение полного материального баланса [7] в ходе итераций и ведёт к сокращению их числа — фактически, снижается размерность пространства, в котором ищется очередное приближение. Ещё одно важное свойство аппроксимации сокиш заключается в том, что его использование гарантирует невырожденность матрицы В предобуславливателя [8]. В свою очередь использование третьего способа (го-№8ит) гарантирует сохранение локального материального баланса. Как показывают результаты численного исследования вариантов сокиш и го'№8ит предобуславливателя МР (раздел 4), эффективность обоих вариантов сравнима.

Для описанного выше примера расчётной сетки 4x4x2 (рис. 1) разложение (1) матрицы предобуславливателя В имеет вид:

В =

Р1 0 0 0

п1 Р 0 0

0 п Р3 0

\ 0 0 пз Р4 /

0

0 0

Р--1Ю1

0 0

0

Р-1Ю2 0

0 0

Р..-1

\

3 ю3

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

Вычисление = В. 1 х происходит аналогично обычному предобуславливанию типа неполного ЬИ разложения и сводится к «прямому ходу» с матрицей ( Р + п)

и затем «обратному ходу» с матрицей ( I + Р 1ю): у находится из уравнения В у = х последовательным решением двух треугольных систем (Р + п)х = х, (I + Р-1ю)у = см. (1).

Таким образом, задача вычисления у = В 1х сводится к вычислению Pi

W

для каждого слоя ячеек. При этом пространственная размерность задачи уменьшается с 3 до 2.

В соответствии с (2) Р^ также представляется в блочном виде Р^ = (Т + т) ■ (1 + Т-1т):

P

( Ti 0 0

т1 T 0

0 т3 T3

\ 0 0 т3

0 0 0

Т4

/

0

0 0

т-У

0 0

0

T-1v2 0

0 0

Т-1

\

3 ^з

где Т^ — матрицы, соответствующие связям в пределах одной «линии» ячеек, слоя с номером %, тг и щ — матрицы, элементы которых соответствуют связям между «линиями» ячеек в пределах слоя .

Процедура вычисления Р-1х полностью аналогична вычислению В-1х. Следующий уровень рекурсии требует обращения трёхдиагональной матрицы Т, которое легко может быть выполнено методом прогонки.

3. Предобуславливатель Generalized Nested Factorization

Предобуславливатель Generalized Nested Factorization впервые рассмотрен в работе [6]. Он является обобщением метода Nested Factorization на случай неструктурированных сеток и наличия нелокальных связей между ячейками.

Основным отличием метода GNF от алгоритма NF является отказ от требования строгой трёхдиагональной структуры (блочной трёхдиагональной структуры для многофазных задач) на самом нижнем уровне построения предобуславлива-теля. Вместо этого предполагается, что центральная блочная полоса матрицы (в оригинальном NF элемент этой полосы — диагональный блок Tj, отвечающий за взаимодействие ячеек внутри одной «линии» ячеек), является плотной матрицей, которая обращается точно. Такие блоки называются «суперклетками». Появившаяся свобода позволяет создавать слои суперклеток, каждая из которых объединяет в себе произвольное число «линий» ячеек по всем трём направлениям пространства. Ячейки должны объединяться в суперклетки таким образом, чтобы, как и в оригинальном NF, матрица предобуславливателя имела на каждом шаге построения блочную трёхдиагональную структуру. В этом случае любой слой (суперклетка) граничит с не более чем двумя соседними слоями (суперклетками), что гарантирует блочную трёхдиагональную структуру на верхнем уровне. Изменение способа объединения ячеек в суперклетки не меняет решения системы, но, в случае GNF, может существенно поменять качество предобуславливания.

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

Очевидно, что, чем больше размер каждой суперклетки, тем более точное решение может быть получено. В предельном случае одной суперклетки (соответствующей всей сеточной расчётной области) предобуславливатель в точности совпадёт с обратной матрицей системы.

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

Здесь каждая суперклетка первого слоя сетки и одна суперклетка второго слоя содержат по две «линии» ячеек, две суперклетки второго сеточного слоя — по одной «линии» ячеек сетки.

Ti V2 W-2

т2 Т2

п2 Тз v4 w3

7П4 Т4 V5

ш5 Т5

Пг Т6 v7

га 7 т7 v&

т8 Т8 V9

т9 Т9

Рис. 2. Сетка (слева) и структура матрицы В для сетки 4 х 4 х 2 при использовании предобуславливателя GNF

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

4. Сравнительный анализ эффективности алгоритма GNF и ILU-предобуславливателей

Тест № 1 — модель «SPE9» Американского общества инженеров нефтяников [9], тест № 2 описан в [10]. Тесты № 3 и № 4 соответствуют реальным месторождениям Западной Сибири.

Используемые матрицы были выгружены из гидродинамического симулято-ра многофазной многокомпонентной фильтрации. Все матрицы были выгружены приблизительно с середины полного времени расчёта модели со второй - третьей нелинейной итерации метода Ньютона. Используемая модель фильтрации — трёхфазная трёхкомпонентная модель «чёрная нефть» [7].

Далее рассмотрим результаты сравнительного анализа эффективности предобуславливателя GNF и предобуславливателей семейства ILU: ILU(0), ILU(1), ILUT — универсальных и широко распространённых в настоящее время [1].

В качестве метода решения СЛАУ был выбран алгоритм FGMRes(k) — «гибкий» вариант метода GMRes(k), позволяющий использовать на каждой итерации алгоритма свой предобуславливатель [11].

Используемые для расчёта матрицы были получены из полных матриц расчётных моделей, описанных выше, с исключением уравнений для скважин с помощью метода дополнения по Шуру [12].

Итерации метода FGMRes(k) продолжались до достижения неравенства ||г||2 ^ £г||г0||2 + ^а при sr = 10-6, £а = 10-8 (|И|2 и ||г0||2 — текущая и начальная невязки, соответственно). Максимальное число итераций равнялось 500. Размерность к крыловского подпространства равнялась 30. В качестве начального приближения выбирался вектор х с компонентами Xi = г.

Параметры предобуславливателя ILUT(/,s) (в строках матриц L и U сохранялось не более I коэффициентов, превышающих по модулю е) имели значения I = 60 и б = 10-4.

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

Обращение внутренних матричных блоков для суперклеток производилось приближённо, с помощью предобуславливателя ILUT (60,10-3).

Результаты исследования времени работы итерационного метода FGMRes с различными предобуславливателями для рассматриваемых моделей представлены в табл. 1.

Таблица 1

Время работы итерационного метода FGMRes(k) с различными предобуславливателями на моделях 1—4. Время указано в секундах, Nu — число итераций, Тргес — время построения предобуславливателя, Тцег — время, затрачиваемое только на итерации, Тац — полное время решения

СЛАУ

Модель № 1 Nlt (с) Tprec (%) Titer (%) Tall (с) Tilut/Tall

ILU(0) 148 32 68 0.38 0.8

ILU(1) 55 51 49 0.31 1.0

ILUT 43 91 9 0.31 1.0

GNF rowsum 15 72 28 0.14 2.2

GNF colsum 16 70 30 0.15 2.1

Модель № 2 Nu (с) Tprec (%) Titer (%) Tall (с) Tilut/Tall

GNF rowsum 142 42 58 0.19 -

GNF colsum 24 50 50 0.12 -

Модель № 3 Nu (с) Tprec (%) Titer (%) Tall (с) Tilut/Tall

ILUT 287 35 65 1.23 1.0

GNF rowsum 30 59 41 0.81 1.5

GNF colsum 38 53 47 0.94 1.3

Модель № 4 Nu (с) Tprec (%) Titer (%) Tall (с) Tilut/Tall

ILU(1) 22 39 61 0.31 0.9

ILUT 17 63 37 0.29 1.0

GNF rowsum 10 46 54 0.15 1.9

GNF colsum 11 41 59 0.16 1.8

Краткое описание исследуемых для каждой модели матриц приведено в таблице 2.

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

В последней графе каждой таблицы приведено отношение «эталонного» времени решения СЛАУ ко времени решения СЛАУ с используемым предобуславли-вателем.

Как видно из табл. 1, предобуславливание методом СМР в целом является эффективным, до 2 раз опережая на всех моделях предобуславливатели семейства 1Ьи по времени работы, и до 9 раз по количеству итераций. Кроме того, для модели № 2 только применение предобуславливателя СМР гарантирует сходимость итерационного метода за максимально отведённое количество итераций.

Таблица 2

Параметры матриц моделей. Обозначения фаз: «н» — нефтяная, «в» —

водная, «г» — газовая

Модель 1 2 3 4

Сетка 24 х 25 х 15 40 х 40 х 1 51 х 51 х 15 22 х 69 х 21

Число активных 9 000 1 600 21 489 31 762

ячеек

Число строк мат- 27 000 4 800 42 978 63 524

рицы

Число ненулевых 542 900 36 043 551 980 423 043

элементов

Число фаз 3 (н,в,г) 3 (н,в,г) 2 (н,в) 2 (в,г)

Число скважин 26 1 2 115

Сравнивая влияние выбора корректировки матрицы ошибки в предобуславлива-теле GNF, можно сделать вывод, что корректировка rowsum обеспечивает в целом лучшую робастность, чем colsum, однако разница во времени работы не является существенной.

В итоге можно сделать следующие выводы:

— «Простые» предобуславливатели ILU(0) и ILU(1) не эффективны для задач рассматриваемого типа.

— Сравнительно хорошие результаты (по числу итераций и времени работы алгоритма) даёт лишь применение предобуславливателя ILUT с достаточно малым порогом отбрасывания ненулевых элементов (ILUT(60,10-4)).

— Наиболее эффективно при использовании рассмотренных тестов ведёт себя предобуславливатель GNF — как по числу итераций, так и по времени работы алгоритма. В случае рассмотренных примеров это связано как с меньшим числом обусловленности предобусловленной матрицы, так и со свойством рассматриваемого предобуславливателя лучше кластеризовать собственные значения (см. раздел 5).

Эффективность использования предобуславливателей NF и GNF, видимо, обусловлена следующими факторами. Скорость построения и применения предобу-славливателя в значительной степени обусловлена рекурсивной трёхдиагональ-ной структурой матрицы В. При этом, если расчётная сетка незначительно отличается от структурированной (например, число неактивных ячеек сетки мало), то и эффективность GNF будет близка к эффективности стандартного NF.

Перенумерация ячеек сетки в соответствии с выбором «главной» координатной оси (см. раздел 2) концентрирует элементы матрицы с большими абсолютными значениями вблизи главной диагонали. Именно эти элементы входят в диагональные блоки Ti матрицы В, которые обращаются точно. Оставшиеся элементы матрицы сравнительно малы. Это повышает качество аппроксимации матрицы А матрицей В. Переупорядочение осей, таким образом, является аналогом выбору ведущего элемента при использовании стандартных предобуславливателей типа ILU — но происходит не поэлементно, а сразу — достаточно большими блоками, соответствующими слоям ячеек сетки. При наличии в задаче нелокальных связей между ячейками (например, скважин) разбиение должно проводиться таким образом, что бы все ячейки, связанные одной скважиной, попадали бы в одну суперклетку. Такая суперклетка, в общем случае, не будет разреженной диагональной матрицей, однако это не является проблемой в силу того, что матрица, соответствующая суперклетке, обращается «целиком». При этом могут быть использованы эффективные прямые алгоритмы для ленточных матриц, например, из пакета LAPACK [13].

Помимо этого, предобуславливатель GNF как в варианте colsum, так и в варианте rowsum, сохраняют важные свойства матрицы исходной системы, связанные

с физикой задачи (см. раздел 2), что также положительно сказывается на скорости сходимости метода.

В отличие от предобуславливателей семейства МР, общие предобуславливате-ли семейства 1Ьи являются «универсальными» в том смысле, что не используют какой-либо априорной информации о структуре математической модели, расчётной сетки и используемых аппроксимациях. Этим объясняется их меньшая эффективность для решения рассматриваемого класса задач.

5. Спектральный анализ предобусловленных матриц

В данном разделе представлены результаты сравнительного анализа спектральных свойств предобуславливателя СМР и предобуславливателей семейства 1Ьи (1Ш(0), 1Ш(1), 1ШТ).

Для вычисления спектров использовался пакет АБРАСК [14]. Реализованные в пакете методы являются методами на основе подпространств Крылова и не требуют, в явном виде, матрицы А или В-1 А.

Известно [1], что скорость сходимости итерационных методов для решения СЛАУ напрямую зависит от числа обусловленности матрицы системы, которое обычно вычисляется в спектральной норме [15]. Для числа обусловленности несимметричной неопределённой системы справедлива оценка

СОП(1(А) > |Атах(А)/Ат;п(А)| ,

где Атт,тах ( А) — минимальное и максимальное собственные числа матрицы А, соответственно.

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

|А| = 1 [15].

Спектры исследуемых предобусловленных матриц представлены на рис. 3-6.

GNF-colsum

GNF-rowsum

Рис. 3. Спектры предобусловленных матриц для модели № 1

GNF-colsum

GNF-rowsum

Рис. 4. Спектры предобусловленных матриц для модели № 2

GNF-colsum

GNF-rowsum

Рис. 5. Спектры предобусловленных матриц для модели № 3

На рисунках показано распределение 500 минимальных и 500 максимальных по модулю собственных значений.

Для моделей № 2 и № 4 результат применения предобуславливателей 1Ьи(0) и 1Ьи(1) не показан, так как в этом случае итерации не сходились.

Полученные результаты исследования спектров хорошо согласуются с численными результатами исследования работы итерационного метода FGMRes(k), во многом объясняя причины численной эффективности алгоритма предобуславли-вания GNF: он обеспечивает лучшую по сравнению с остальными предобуслав-ливателями кластеризацию спектра и число обусловленности предобусловленной матрицы системы, что наиболее заметно для модели № 3.

GNF-colsum GNF-rowsum

Рис. 6. Спектры предобусловленных матриц для модели № 4

6. Заключение

Данная работа посвящена исследованию метода предобуславливания Generalized Nested Factorization (GNF) для реальных задач пластовой фильтрации.

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

По сравнению с универсальными ILU - предобуславливателями, являющимися на сегодняшний день стандартными для решаемого класса задач, время решения одной СЛАУ уменьшается до 2 раз, количество итераций сокращается до 9 раз, что позволяет сделать вывод о хорошей робастности предобуславливателя GNF для реальных задач пластовой фильтрации.

Литература

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

2. Баландин М. Ю., Шурина Э. П. Методы решения СЛАУ большой размерности. — Новосибирск: Издательство НГТУ, 2000. [Balandin M. Yu, Shurina A. P. Solution of Large Linear Systems — Novosibirsk: Publisher NSTU, 2000 ]

3. Benzi M. Preconditioning Techniques for Large Linear Systems: A Survey // Journal of Computational Physics. — 2002. — Vol. 182, No 2. — Pp. 418-477.

4. Constrained Residual Acceleration of Conjugate Residual Methods / J. R. Wallis, R. P. Kendall, T. E. Little, J. S. Nolen // SPE Paper 13536, 8th SPE Symposium on Reservoir Simulation, Dallas, TX. — 1985.

5. Appleyard J. R, Cheshire I. M. Nested Factorization // SPE Paper 12264-MS, SPE Reservoir Simulation Symposium. — No 12264, 12264. — San Francisco, California: Society of Petroleum Engineers, 1983. — Pp. 315-324.

6. Appleyard J. R. Method and Apparatus for Estimating the Physical State of a Physical System // European Patent Application EP 2 068 263 A2. — 2010.

7. Aziz K., Settari A. Petroleum Reservoir Simulation. — Essex, England: Elsevier Applied Science Publishers, 1979.

8. Fourier Analysis of Modified Nested Factorization Preconditioner for Three-Dimensional Isotropic Problems: Research report / P. Kumar, L. Grigori, Q. Niu, F. Nataf. — 2010. — http://hal.inria.fr/inria-00448291.

9. Killough J. E. Ninth SPE Comparative Solution Project: A Reexamination of Black-Oil Simulation // SPE Paper 29110-MS, SPE Reservoir Simulation Symposium. — San Antonio, Texas: Society of Petroleum Engineers, 1995. — Pp. 135147.

10. Некрасов А. А. Тесты и методики испытания программного обеспечения моделирования разработки месторождений // Вестник ЦКР Роснедра. — 2008. — № 1. — С. 54-64. [Nekrasov A. A. Testing Methods of Simulation Software for Field Development // Vestnik CKR Rosnedra. — 2008. — No. 1. — P. 54-64 ]

11. Saad Y. A Flexible Inner-Outer Preconditioned GMRES Algorithm // SIAM Journal on Scientific Computing. — 1993. — Vol. 14, No 2. — Pp. 461-469.

12. Прасолов В. В. Задачи и теоремы линейной алгебры. — М.: Физматлит, 1996. [Prasolov V. V. Problems and Theorems of Linear Algebra. — Moscow: Fizmatlit, 1996 ]

13. LAPACK Users' Guide (3nd edition) / E. Anderson, Z. Bai, C. Bischof et al. — Philadelphia: SIAM, 1999.

14. Lehoucq R., Sorensen D., Yang C. ARPACK Users' Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. — Computer Based Learning Unit, University of Leeds, 1998.

15. Деммель Д. Вычислительная линейная алгебра. Теория и приложения. — Москва: Мир, 2001. [Demmel J. Applied Numerical Linear Algebra. Theory and Applications. — Moscow: MIR, 2001 ]

UDC 519.6

Numerical Investigation of Generalized Nested Factorization Preconditioner for Reservoir Simulation Problems

V. E. Borisov*, E. B. Savenkov^

* Department of mathematics and mechanics Moscow State University

Leninskiye Gory, 1, Moscow, Russia, 119991 t Keldysh Institute for Applied Mathematics RAS Miusskaya Sq., 4a, Moscow, Russia, 125047

In this paper we present some results on application of Generalized Nested Factorization (GNF) preconditioner for real-life reservoir simulation problems. The test problems considered in the paper share some important features with reservoir models arise in real-life practice, e.g., a presence of non-local cell connections and unstructured computational grids.

Numerical results are compared with the ones obtained using standard ILUT preconditioning techniques. Spectral properties of preconditioned matrices are analyzed.

As a result very good performance of GNF preconditioner is observed for real-life reservoir simulation problems.

Key words and phrases: reservoir simulation, preconditioning, incomplete LU factorization, Krylov subspace methods, generalized nested factorization.