CC BY
748 вестн. моек. ун-та. сер. 1, математика. механика. 2009. №б
Пусть теперь последовательность Л ограничена. Тогда сходимость ряда YlkLi §¡T эквивалентна сходимости ряда £a,. По теореме C получаем, что любая функция, для которой M,[f,I] ^ a, при всех к, имеет ограниченную вариацию, т.е. попадает в класс ЛBV(I) = BV(I).
Необходим остъ. Пусть последовательности {а^} и {Ад,} таковы, что ^fcli §¡; = Обозначим отрезок I через [bo, b] и возьмем на нем строго возрастающую последовательность точек {b,}, стремящуюся к b.
Рассмотрим вначале случай a, — О. Тогда ряд ^j=l(— 1)j-1aj сходится по признаку Лейбница, обозначим его сумму через a. Пусть f (bo) = О, f (b= Ylj=l(- 1)j-1aj, на [b,-l,b,] зададим функцию линейным образом, и пусть f (b) = a. Мы получили непрерывную на отрезке функцию. С одной стороны, рассматривая системы интервалов {(bj-l,bj)^=1, мы видим, что Л-вариация построенной функции бесконечна. С другой стороны, рассматривая функции
gk(x) = ( f (x), x ^ h; [f(b,), x > b,,
которые имеют к промежутков монотонности, мы видим, что sup \ f—g,\ = \ f (b,+l) — g,(b,+l)\ = a,+l ^ a,, т.е. функция f имеет подходящий порядок убывания кусочно-монотонных приближений.
Пусть теперь a, — ao > О. Положим тогда f (b, ) = (—ao, на [b,-l,b, ] зададим функцию линейным образом, и пусть f (b) = О. Тогда M,[f,I] ^ ao ^ a,, что сразу видно, если взять приближающую функцию g,(x) = О, но в то же время f не принадлежит ни одному классу Ватермана, поскольку не имеет предела слева в точке b. Теорема доказана.
Автор признателен профессору М.И. Дьяченко за обсуждение работы.
Работа выполнена при финансовой поддержке РФФИ (проект № 08-01-00302) и программы "Ведущие научные школы" (проект НШ-2787.2008.1).
СПИСОК ЛИТЕРАТУРЫ
1. Waterman D. On convergence of Fourier series of functions of generalized bounded variation // Stud. Math. 1972. 44, N 1. 107-117.
2. Чантурия З.А. Модуль изменения и его применение в теории рядов Фурье // Докл. АН СССР. 1974. 214, № 1. 63-66.
3. Севастьянов Е.А. Кусочно-монотонная аппроксимация и Ф-вариация // Anal. Math. 1975. 1. 141-164.
4. Avdispahic M. On the classes ЛBV and V[v] // Proc. Amer. Math. Soc. 1985. 95, N 2. 230-234.
Поступила в редакцию 28.11.2008
УДК 519.6
ПРИМЕНЕНИЕ ПАРАЛЛЕЛЬНОГО ПРЕДОБУСЛОВЛИВАТЕЛЯ CPR К ЗАДАЧЕ ФИЛЬТРАЦИИ ВЯЗКОЙ СЖИМАЕМОЙ ЖИДКОСТИ
В ПОРИСТОЙ СРЕДЕ
И. Г. Горелов1
Рассмотрен метод предобусловливания СРИ применительно к решению систем с матрицей, возникающих при аппроксимации систем дифференциальных уравнений в частных производных для задач фильтрации многокомпонентной смеси в пористой среде. Предложены новые версии алгоритма. Проведены численные эксперименты в составе российского промышленного комплекса гидродинамического моделирования на параллельных ЭВМ
1 Горелов Илья Геннадьевич — асп. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: chronosphere@mail.ru.
ВЕСТН. МОСК. УН-ТА. СЕР.1, МАТЕМАТИКА. МЕХАНИКА. 2009. №6
49
для тестовых и реальных моделей газонефтяных пластов Западной Сибири, подтверждающие эффективность методов.
Ключевые слова: CPR-метод, параллельный предобусловливатель, разреженные матрицы, итерационный алгоритм, течение вязкой жидкости в пористой среде.
The constrained pressure residual (CPR) preconditioning method is considered with regard to solution of systems with matrices appearing in discretization of PDE systems describing multicomponent fluid flow in porous media. New versions of algorithms are proposed. Numerical experiments using an actual parallel hydrodynamic simulator were performed for imitation and actual oil fields in Western Siberia, which confirm the efficiency of the methods.
Key words: CPR-method, parallel preconditioner, sparse matrices, iteration algorithm, fluid flow in porous media.
Рассмотрим систему уравнений в трехмерной области [1]: д к
— (фNc) = div ^ xCtpCpl3{k^(Vpp - ppgVD)) + qc, P=O,W,G HP
С = 1,..., Пс , Пс = 3, PO - PG = Pcog, PO - PW = Pcow,
Sw
Sw + So + Sg = 1, Ni = Nw = £w,
SC
Bw '
N2 = No = ZoMlT + RO,GN3 = Ng= ZaMir +
bo bg bg bo
Искомыми величинами являются:
Nc = Nc(t, x, y, z) — молярные плотности; pP = pP(t, x, y, z), P = O,W,G — давление; Sp = Sp(t,x,y,z) — насыщенности фаз.
Величины ф, xc,P, £p, k, krP, hp, pP, qc, Pcog, Pcow, Bp, Rg,o, Ro,g, являющиеся заданными функциями, нелинейно зависящими от искомых величин, определены в [1]; в, g и sc — известные константы, а D — известная функция. Для аппроксимации по времени используется полностью неявная схема (см. [1]). Аппроксимацией по пространственным переменным методом конечных объемов [1] задача сводится к системе нелинейных алгебраических уравнений вида F(p, Ni,..., N„c) = 0, где p = (p%), Nc = (NC) — векторы значений давления и молярных плотностей в блоках сетки. Для решения системы используется метод Ньютона. На каждом шаге метода Ньютона необходимо решать систему с несимметричной матрицей — якобианом из метода Ньютона Ax = b, где A можно рассматривать как матрицу, элементами которой являются блоки размера (пс) х (пс).
Отыскание предобусловливателя Mcpr заключается в построении по матрице A подматрицы Ap, при котором используются элементы блоков, соответствующие только переменной давления. Действие Mcpr выглядит так [2]: решаем систему
Apxp — rp, (1) по xp строим вектор xz (вставляя нули на места неизвестных) и подставляем его в
Mcpr r = Mr + xz - MAxz = xcpr.
Здесь M — любой "обычный" предобусловливатель.
Предположим, что вместо элементов, соответствующих переменной давления, мы используем для построения элементы, выбираемые исходя из некоторого иного критерия. Будем называть такой метод Adaptive CPR или ACPR (PACPR в параллельном случае). Множество переменных обозначим через T, а искомую матрицу — через At. За количество элементов, которые будут участвовать в построении,
50
ВЕСТН. моек. УН-ТА. СЕР.1, МАТЕМАТИКА. МЕХАНИКА. 2009. №6
отвечают Ei £ [0%, 100%]: E\ соответствует давлению, а и E3 — молярным плотностям. Отметим, что при значениях параметров (Ei,E2,E3) = (100, 0, 0) мы получаем классический алгоритм CPR.
В качестве критерия отбора элементов возможно использовать, например, принадлежность множеству блоков фронта. Под фронтом при этом понимаются блоки, значения неизвестных в которых меняются наиболее сильно на данном временном шаге. Для приближенного определения таких блоков можно произвести одну или несколько итераций обычного алгоритма решения линейной системы, а потом воспользоваться полученным грубым решением.
Многопроцессорная версия описанного метода была построена с применением перекрытия расчетных областей в разных процессах [3].
Такой подход был реализован в виде параллельного программного модуля и протестирован в составе российского промышленного комплекса гидродинамического моделирования. В качестве тестовых были выбраны модели реальных месторождений Западной Сибири и несколько синтетических (признанных сообществом инженеров-нефтяников как эталонные для тестирования работоспособности гидродинамического симулятора). Численные эксперименты показали, что алгоритм работает быстро, но высокая скорость достигается при разных значениях параметров метода при переходе от модели к модели.
В ходе тестирования было установлено, что для работоспособности CPR достаточно очень грубого приближения к точному решению системы (1). Поэтому может оказаться, что достаточно не решать систему, а просто хорошо предобусловить ее. При таком подходе отпадает необходимость в хранении самой матрицы давления и части дополнительных массивов, а также появляется возможность равномерно распределить нагрузку между вычислительными ядрами и объединить часть операций.
В качестве локального предобуслов-ливателя возьмем метод ILU(0), а глобального — блочный алгоритм того же метода. Шаблон (множество ненулевых элементов) матрицы давления в точности совпадает с таковым для исходной матрицы блоков. Это позволяет нам за один проход по матрице A построить сразу два предобусловливателя, а также не строить и не хранить дополнительный массив с номерами ненулевых элементов. Приведенный вариант CPR назовем BCPR (и соответственно PBCPR). Он также был реализован в виде параллельного программного модуля в составе российского промышленного комплекса гидродинамического моделирования. Результаты тестов на машине с восемью ядрами частотой 2,8 ГГц и общей памятью 16 Гб показаны в таблице. Числитель дроби соответствует методу решения с использованием параллельного блочного неполного LU-разложения, а знаменатель — параллельного блочного метода CPR. Приведено время всего расчета модели, включающее чтение данных, построение якобиана, решение системы и так далее.
Из таблицы видно, что PBCPR работает немного быстрее базового алгоритма блочного параллельного ILU(0) (являющегося стандартом для данного программного пакета). Более детальное тестирование показывает, что множество матриц, на которых алгоритм сохраняет работоспособность, практически совпадает с таковым для исходного блочного параллельного ILU(0). Все методы хорошо масштабируются при увеличении числа расчетных ядер.
Построенный алгоритм выглядит перспективным с точки зрения промышленного применения на рассматриваемом классе задач.
СПИСОК ЛИТЕРАТУРЫ
1. Aziz K., Settari A. Petroleum Reservoir Simulation. London: Applied Science Publishers, 1979.
2. Wallis J. R, Kendall R. P., Little T. E. Constraint Residual Acceleration of Conjugate Residual Methods // Reservoir Simulation Symposium. Dallas, Texas, February 10-13, 1985. SPE 13536.
3. Капорин И.Е., Коньшин И.Н. Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки // Журн. вычисл. матем. и матем. физ. 2001. 41, № 4. 515-528.
Поступила в редакцию 29.04.2009
Модель Число Время, с
неизвестных итераций метода BCGStab итераций метода Ньютона
1 (486 х 3) 4224/3378 356/305 1/1
2 (31923 х 3) 485/292 25/25 2/2
3 (9000 х 3) 1994/1476 169/168 8/7
4 (255366 х 3) 737/437 102/102 124/115
5 (255366 х 3) 714/347 76/76 162/153
6 (2307383 х 3) 1739/1150 304/304 3939/3715
