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

На правах рукописи

4852017

ГОРЕЛОВ Илья Геннадьевич

Применение параллельного

предобусловливателя СРЫ к задаче фильтрации вязкой сжимаемой жидкости в пористой среде

05.13.18 — математическое моделирование, численные методы и комплексы программ

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

КАЗАНЬ 2011

(1 1 АВГ 2011

4852017

Работа выполнена на кафедре вычислительной математики механико-математического факультета

Московского государственного университета им. М.В.Ломоносова

Научный руководитель:

кандидат физико-математических наук, доцент Богачев Кирилл Юрьевич.

Официальные оппоненты: доктор физико-математических наук,

профессор Лапин Александр Васильевич;

кандидат физико-математических наук, доцент Вершинин Анатолий Викторович.

Ведущая организация:

Институт прикладной математики им. М.В. Келдыша РАН

Защита диссертации состоится 22 сентября 2011 г. в 14 часов 30 минут на заседании диссертационного совета Д 212.081.21 при ФГАОУВПО «Казанский (Приволжский) федеральный университет» по адресу: 420008, г. Казань, ул. Кремлевская, д.18, корп.2, ауд. 218.

С диссертацией можно ознакомиться в научной библиотеке им. Н. И. Лобачевского ФГАОУВПО «Казанский (Приволжский) федеральный университет».

Автореферат разослан «19» августа 2011 г.

Ученый секретарь диссертационного совета, доктор физико-математических наук,

профессор

Задворнов О.А.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

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

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

Существует множество способов решения систем линейных алгебраических уравнений с разреженной матрицей и их предобусловливания, рассмотренных такими авторами как Г. И. Марчук, Н.С. Бахвалов, Г.М. Кобельков, С.К. Годунов, Р.П. Федоренко, И. Саад и другими. Однако лишь очень немногие из них учитывают специфику данной конкретной задачи (например методы, описанные в работах Д. Уолиса, К. Азиза, К. Штубена). Между тем, именно специфические алгоритмы, использующие различные особенности исходной системы дифференциальных уравнений, практически всегда оказываются наиболее эффективными.

Одним из таких алгоритмов является метод CPR (Constrained Pressure Residual), основанный на идее использования при построении предобусловливателя не только исходной системы линейных алгебраических уравнений, но и матрицы, фактически отвечающей за неявную только по давлению аппроксимацию и потому имеющей в несколько раз меньшую размерность (в 3 раза для трехкомпонентной смеси). Существует множество как академических, так и прикладных публикаций, в которых утверждается, что этот метод, при его использовании совместно с алгебраическим многосеточным алгоритмом, является более эффективным по сравнению с предобусловливателями «общего назначения» (например, неполное LU-разложение или ILU). Причем он позволяет ускорить расчет не на несколько процентов, а в несколько раз. Тем не менее, на данный момент практически ни

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

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

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

при промышленных расчетах.

Цели работы и задачи исследования

1. Анализ алгоритма СРИ применительно к решению систем линейных алгебраических уравнений, возникающих при моделировании процессов фильтрации вязкой сжимаемой жидкости в пористой среде с учетом факторов, характерных для реальных задача и оказывающих сильное влияние на получающуюся в результате матрицу и, как следствие, па сходимость метода.

2. Разработка на его базе новых, более эффективных с точки зрения скорости работы и универсальности применения алгоритмов.

3. Разработка параллельных программных модулей, реализующих методы.

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

5. Определение оптимальных параметров алгоритмов и степени влияния того или иного свойства пласта и гидродинамической модели на сходимость.

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

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

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

Научная новизна

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

2. Предложены новые, в том числе более эффективные, версии алгоритма, учитывающие, в частности, архитектуру современных

ЭВМ.

3. Разработаны программные модули, реализующие параллельные версии всех рассмотренных методов, которые показали свою применимость при решении реальных задач.

4. Определены преимущества и недостатки алгоритмов, подтвержденные экспериментальным путем при проведении численного моделирования на параллельных ЭВМ с использованием тестовых и реальных гидродинамических моделей газонефтяных пластов. Подтверждена эффективность методов и описан круг задач, для которых использование того или иного способа решения является наиболее оптимальным.

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

Положения, выносимые на защиту

1. На основе анализа метода СРИ предложены более эффективные и обобщающие версии алгоритма, пригодные для использования при

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

2. Разработан подход к программированию, учитывающий особенности предлагаемых методов и архитектуру современных многопроцессорных ЭВМ и позволяющий добиться эффективной работы.

3. Разработаны программные модули, реализующие параллельные версии рассмотренных методов.

4. Выявлены теоретически и подтверждены экспериментально путем численного моделирования на параллельных ЭВМ с использованием тестовых и реальных гидродинамических моделей газонефтяных пластов различные свойства и параметры алгоритмов. Описал круг задач, для которых применение этих методов наиболее оправдано.

Апробация работы. Основные результаты докладывались и обсуждались на:

• научных семинарах каф. вычислительной математики механико-математического факультета МГУ им. М.В. Ломоносова, Москва, 2007-2009 гг.,

• научном семинаре Института вычислительной математики РАН, Москва, 2008 г.,

• научном семинаре кафедры дифференциальных уравнений Московского энергетического института, Москва, 2009 г.,

• конференциях «Ломоносовские чтения», Москва, 2008, 2009 гг.,

• 4-й международной конференции «Математические идеи П.Л. Чебышева и их приложения к современным проблемам естествознания», Обнинск, 2008 г.,

• научном семинаре кафедры математического моделирования МГТУ им. Н.Э.Баумана, Москва, 2009 г.

• научном семинаре кафедры вычислительной математики Казанского (Приволжского) федерального университета, Казань, 2011 г.

Публикации. Научные результаты опубликованы в 4 работах, в том числе трех статьях в журналах перечня ведущих рецензируемых научных журналов и изданий ВАК РФ [1, 2, 3]. Журнал [4] вошел в перечень ВАК уже после выхода публикации.

Личный вклад соискателя. Все исследования, изложенные в диссертационной работе и выносимые на защиту, проведены лично соискателем. Во время исследований были использованы визуализаторы, загрузчики и другие функциональные блоки, являющиеся плодом коллективного труда студентов, аспирантов и преподавателей Московского Государственного Университета им. М.В. Ломоносова.

Структура и объем работы. Диссертационная работа состоит из введения, четырех глав, выводов, заключения, приложения и списка литературы. Полный объем составляет 104 страницы. Библиография включает 29 наименований.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ

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

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

Исходная задача сводится к следующей системе уравнений (стандартная трехфазная модель черной нефти, Коновалов, Азиз):

+ Яс,

с=1 ,...,пс, пс = 3,

РО - Рсг = Рсов, РО~Р№ = РС0\У, + + 5с = 1,

Ач = А» = Ы.зс^г,

(1)

Искомые величины:

1. Nc = Nc(t,x,y,z) — молярная плотность компонента с. Для модели черной нефти компонентами служат вода, условная нефть (тяжелые углеводороды) и условный газ (легкие углеводороды), для композиционной модели — различные составные фракции перечисленных веществ;

2. Pw = Pw(t,x,y,z),p0=p0(t,x,y,z),pG = PG(t,x,y,z) - давление водяной, нефтяной, газовой фаз соответственно;

3. Sp = Sp(t,x,y,z) — насыщенность Р-ой фазы, Р = 0,G, W, — выражается явно через давление и молярные плотности с помощью последних формул системы.

Обозначения:

• Р = О, G, W — фазы: нефть, газ, вода;

• р — давление в одпой из фаз: в нефтяной фазе (ро), если в модели есть нефть, иначе в газовой (pg);

• N = (Ni,..., Nnc) вектор молярных плотностей;

• Rg,o — Rg,o{p> N) — растворимость газа в нефтяной фазе;

• Ro,g ~ Ro.ciP, N) — летучесть нефти в газовой фазе;

• Bp — Вр(р, N) — коэфф. объемного расширения фазы;

• ф = ф(р, х, у, z) — пористость среды;

• £с,р = £Cip(p, N) — молярная доля компонента с в фазе Р;

• £р = £р(р> N) — молярная плотность фазы;

• к = к(р, х, у, z) — тензор абсолютной проницаемости;

• kr%p = kTip(Sw, SG) — относит, фазовая проницаемость;

• fip = Цр(р, N) — вязкость фазы;

• D = D(x,y, z) — глубина (сверху вниз);

• рр = рр(р, N) — массовая плотность фазы;

• Рсоа = Pcog{Sg), Pcow = Pcow{Sw) — капиллярное давление в системе нефть-газ, в системе вода-нефть;

9 Чс — Чс(р, t, х, у, z, N) — источник компонента с (например скважина);

• д = const — значение ускорения свободного падения;

• (3 = const — константа Дарси;

• Zp,sc = const — молярн. плотн. фаз в поверхностных усл.

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

где р = (р1), 1ЧС (Дгс') векторы значений давления и молярных плотностей в блоках сетки.

Для решения системы нелинейных уравнений ^(С/) = 0, и = (р, ЛГ) используется метод Ньютона:

Здесь 8F(Um)JdU - отображение (матрица) Rn^K'+J) RпМк+-П х

каждом шаге метода Ньютона надо решать систему с несимметричной

F(p,Nb...Nj=0,

Rncx(/f+j), К число блоков сетки, J число скважин. На

матрицей dF(Um)/dU, т.е. задача сводится к решению системы линейных алгебраических уравнений:

Ах = г, (2)

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

Полученная система решается методом BCGStab (Van der Vorst, 1992) с применением предобусловливания с помощью неполного LU разложения. Развитием данной совокупности алгоритмов является внедрение метода CPR. Суть его заключается в использовании особым образом при предобусловливании элементов матрицы А, соответствующих переменным, отвечающим за давление в блоках расчетной сетки.

Преобразованиями матрицы (например, методом Гаусса) обращаем в ноль элементы в диагональных блоках. Для матрицы из четырех блоков

/ au ®15 «iß ^ ^

размером 3x3 это выглядит следующим образом:

f l an аг2 an \ a2i а22 Ö23

у о31 а32 азз }

I

«41 ^42 «43 ^

151 Я52 Ö53

\ «61 «62 Q63 /

024 «25 <126

V «34 «35 азб /

/ Й41 Я45 а.1С ^

054 Й55 «56

У а64 а65 абб J

( а",

п О О

a'2i а'22 О

\ «з 1 «32 «зз /

/ л" Л» \

«41 «42 а 43 \

«и а'52 «'53

V \ «61 «62 «63 /

Далее строим матрицу Ар:

/ л" \ \

а и а 15 а 16 \

г.' г,1

а 24 «25 «26

\ Ö34 «35 «36 /

1 а" 44 О 0 \

а'54 а'55 О

У «64 «65 «CG У

/

а"п О О

\

а'21 а'22 О

V аз1 «32 «зз /

/ л" Л» \

а 41 «42 « 43 I

«'51 «'52 «'оз

14 а 15 « 16 «'24 «'25 «'26 \ а34 а35 а3б )

а" 44 О 0 \ «'54 «'55 О

а и «и «"41 «"44

— Av

\ \ «61 «62 «63 / \ «64 «65 «66 ) ) Эту матрицу можно считать приближением для матрицы, получающейся в схеме, неявной по давлению и явной по остальным переменным.

Действие предобусловливателя Мфг выглядит следующим образом (Дурловский, Азиз):

М^г = Мг + хг- MAxz,

АрХр = Tpj

где M - любой «обычный» предобусловливатель (например ILU(O)), под гр подразумевается сужение, аналогичное используемому при построении Ар, а под .т. результат обратного преобразования.

Модификация метода CPR для параллельных ЭВМ базируется на перекрытии расчетных областей в разных процессах (Капорин, Коныпин). Обозначим

• Xi = (pi, Ni, ■ ■ ■ , i-я компонента решения,

• Хк = - ,xtk^) компоненты, отнесенные к к-иу процессу (штоку исполнения),

• Ак, Ьк — соответствующие части матрицы А (размера тпк х пс х (ff + J)) и правой части Ь в fc-м процессе (потоке).

s

Для алгоритмов с перекрытием Ylmk > (К + J), где s — число

fc=i

процессов. Метод CPR применяется к матрице Ак в каждом процессе:

Mcpr,krk = Mkrk + £к _ мкАкхк = (3)

Акхк = гк, (4)

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

Заметим, что решение системы (4) происходит независимо в каждом из процессов и не требует обмена данными между ними.

Глава 2. Вторая глава посвящена развитию метода CPR. В начале приведено краткое обоснование теоретической возможности его использования совместно с алгебраическим многосеточным методом. Кроме того, описаны теоретические предпосылки возможной нестабильной работы алгоритма алгебраического мяогосеточного

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

Во второй части главы предлагаются два варианта метода CPR.

Первый из них основывается на замене подматрицы, строящейся по неременным, соответствующим давлению, матрицей, строящейся по переменным, выбираемым исходя из некоторого иного критерия. В частности, из критерия принадлежности блока, связанного с переменной, множеству (назовем его Г) блоков так называемого фронта. Описан способ выбора блоков фронта. Алгоритм назван Adaptive CPR или ACPR и строится таким образом, что при некоторых его параметрах мы получаем исходный CPR, а при других — варианты тандема BS-GStab+ILU. Среди основных преимуществ предлагаемого метода можно выделить следующие:

• Универсальность — в предельных случаях получаются известные методы.

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

• Сравнительно простая процедура интеграции в уже имеющийся алгоритм расчета.

• Проверенная на практике хорошая сходимость и устойчивость к нефизичным входным данным в предельных случаях.

• Работа исключительно на уровне матрицы — нет необходимости в какой-либо дополнительной информации о структуре сетки.

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

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

Несмотря на вышеописанные преимущества методы СР11 и АСРИ. страдают от ряда фундаментальных недостатков, делающих их, как показывают численные эксперименты, весьма узкоспециализированными. Среди основных можно выделить следующие:

• Необходимость отдельного построения двух независимых предобус-ловливателей и хранения их и матрицы давления в памяти ЭВМ. Как следствие, большой расход оперативной памяти и замедление алгоритма.

• Необходимость решения системы с локальной матрицей при каждом применении глобального предобусловливателя.

• Наличие параметров, требующих подбора вручную. Более того, как показали проведенные расчеты, не существует «универсальных» значений параметров алгоритмов, оптимальных для всех матриц.

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

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

Как показывают численные эксперименты, для корректной работы алгоритма СРП достаточно выбрать очень низкую точность решения системы (3). Поэтому вместо решения системы может оказаться достаточным просто хорошо иредобусловить ее. Это избавит нас от необходимости подбирать локальную точность. В случае, если в качестве предобусловливателя воспользоваться методом 1И1(0), то, кроме того, мы Еге обретем повых параметров, которые присущи многим другим алгоритмам. Важно подчеркнуть, что такой подход очень выгоден еще и тем, что у нас отпадает необходимость в хранении самой матрицы давления Ар и большинства сопутствующих решению локальной системы массивов — достаточно лишь построенного по ней предобусловливателя. Как следствие, расход памяти существенно сокращается. Кроме того, появляется возможность для ускорения алгоритма объединить многие операции, что позволяет эффективнее использовать локальность данных и повысить загрузку центрального процессора. Еще одним способом ускорить расчет является использование блочного алгоритма 1Ы1(0) в качестве глобального предобусловливателя и обычного — в качестве локального. Основным преимуществом данного подхода является то, что шаблон (множество ненулевых элементов) матрицы давления в точности совпадает с таковым для исходной матрицы блоков. Поэтому и структуры получаемых 1Ш-разложений также будут совпадать. Это позволяет построить 1Ьи(0)-разложение локальной матрицы одновременно с построением такового для матрицы всей системы, что приводит к экономии на обращениях в память, а также более эффективному использованию процессорного кэша. Кроме того, отпадает необходимость строить и хранить дополнительный массив с номерами ненулевых элементов. При этом само по себе использование блочного алгоритма

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

Глава 3. Третья глава посвящена особенностям практической реализации алгоритмов. Здесь изложены используемые методы хранения матриц. В частности, описаны преимущества формата MSR (Modified Sparse Row) над CSR (Compressed Sparse Row) для данной задачи. Рассмотрены особенности параллельных реализаций, способы наиболее оптимального выделения оперативной памяти для их работы и грамотного использования других ресурсов ЭВМ. В завершении главы поэтапно и подробно описана процедура решения системы (2) с использованием параллельной ЭВМ для каждого из рассмотренных алгоритмов. При этом детально изложен ряд способов оптимизации, позволяющих ускорить расчет.

Глава 4. Четвертая глава посвящена результатам численного моделирования процессов фильтрации вязкой сжимаемой жидкости в пористой среде с использованием неструктурированной сетки и учетом неоднородностей среды, наличия внешних источников и стоков, капиллярных сил, гравитации и других факторов, оказывающих негативное влияние на сходимость. Тестирование проводилось с использованием десяти различных гидродинамических моделей. Среди них присутствуют как синтетические тесты, так и модели реальных месторождений Восточной и Западной Сибири. При этом число неизвестных в получающейся системе линейных уравнений для некоторых из них достигает 7 х 106 штук, скважин — нескольких тысяч, а промежуток времени, для которого производится моделирование, равен нескольким годам. Приведены рисунки, на которых изображены

расчетные сетки и отмечено расположение скважин, описаны особенности каждой из моделей.

Далее проанализированы результаты тестирования всех рассмотренных алгоритмов, в том числе ОРИ. вместе с алгебраическим миогосеточным методом. Рассмотрен вопрос о поиске оптимальных значений параметров для СР11 и АСРН. Кроме того, рассмотрены случаи предельных значений параметров для ЛСРИ, проведено исследование влияния на суммарную скорость расчета и сходимость таких факторов, как число производимых итераций для поиска множества блоков фронта и наличие/отсутствие в модели скважин. Приведены результаты тестов масштабируемости параллельных алгоритмов СРП, АСРЯ и ВСР11.

Тестовый стенд, на котором производились все вычисления, представлял из себя машину с двумя 4-х ядерными процессорами частотой 2.8 ГГц с общей памятью объемом 16 Гбайт.

Приложение. В приложении содержатся графики с результатами расчетов.

ОСНОВНЫЕ ВЫВОДЫ И РЕЗУЛЬТАТЫ РАБОТЫ

Параллельная версия алгоритма АСР11 (и ее частный случай — СРЛ) демонстрирует свою эффективность на большом числе реальных промышленных задач. Однако, для его применения необходимо знать набор оптимальных параметров, который, как показывают численные эксперименты, меняется при переходе от одной гидродинамической модели к другой. Данный метод может быть полезен в том случае, когда требуется многократно производить одно и то же моделирование, внося небольшие изменения в модель (добавление скважины, изменение

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

Что касается параллельного метода ВСРИ, при правильной реализации с использованием описанных в диссертационной работе оптимизаций он обладает следующими основными преимуществами:

• в отличие от параллельных СРЛ и АСРН не имеет регулируемых параметров;

• работает быстрее базового блочного ВС051аЬ+В1Ы](0) и других рассмотренных алгоритмов;

• стабильность работы такая же, как и у базового блочного алгоритма;

• уменьшение точек синхронизации (в сравнении с базовым) и более равномерное распределение нагрузки между потоками (по сравнению с СРИ и АСРЛ) положительно сказываются на общей сходимости метода и его масштабируемости;

• размер требуемой памяти незначительно превосходит таковой для блочного метода неполного Ш разложения.

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

Таким образом параллельный алгоритм ВСРЛ по совокупности параметров превосходит все другие рассмотренные методы и может с успехом применяться для гидродинамического моделирования в

программах для промышленных расчетов. Просуммируем основные результаты.

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

2. Разработаны реализации алгоритма СРК, адаптированные под архитектуру современных ЭВМ и позволяющие добиться эффективной работы.

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

4. Проведено моделирование на специальных тестовых задачах, показавшее совпадение результатов с известным решением.

5. По результатам проведенного численного моделирования па параллельных ЭВМ с использованием тестовых и реальных гидродинамических моделей газонефтяных пластов определены преимущества и недостатки методов. Описан круг задач, для которых применение того или иного способа решения является наиболее оптимальным.

СПИСОК РАБОТ АВТОРА ПО ТЕМЕ ДИССЕРТАЦИИ

[1] Горелов И.Г. Применение параллельного предобусловливателя СРИ к задаче фильтрации вязкой сжимаемой жидкости в пористой среде // Вестник МГУ, Сер. 6, Математика, Механика, 2009.

[2] Богачев К.Ю., Михалева М.Ю., Горелов И.Г. Применение алгебраического многоуровневого метода АМС: сравнение с методом ВГССБТАВ + 1Ы] и использование АМС в составе метода СРИ. // Вестник МГУ, Сер. 4, Математика, Механика, 2010.

[3] Горелов И. Г. Параллельный адаптивный предобусловливатель СРП в задаче фильтрации // Программные продукты и системы, номер 3, 2009 г; стр. 107-109.

[4] Богачев К.Ю., Горелов И. Г. Применение параллельного предобусловливателя СР11 к задаче фильтрации вязкой сжимаемой жидкости в пористой среде // Вычислительные методы и программирование, 2008, том 9, стр. 184-190.

Подписано в печать:

18.07.2011

Заказ № 5744 Тираж -100 экз. Печать трафаретная. Типография «11-й ФОРМАТ» ИНН 7726330900 115230, Москва, Варшавское ш., 36 (499) 788-78-56 www.autoreferat.ru

Введение

1 Постановка задачи

1.1 Исходные формулы.

1.2 Аппроксимация.

1.3 Сведение к системе и имеющийся метод решения.

1.3.1 Сведение к системе.

1.3.2 Метод «¡Ьи+ВССБЪаЬ».

1.4 Классический метод ОРЯ.

1.4.1 Идея метода СРЫ

1.4.2 Алгоритм метода СРЯ.

1.4.3 Построение алгоритма СРЫ для параллельных ЭВМ (РСРЯ)

2 Расширения метода СРЯ,

2.1 Применение СРЯ вместе с алгебраическим многосеточным алгоритмом (АА/Ю).

2.2 Адаптивный СРЯ (АСРЯ).

2.2.1 Алгоритм метода АСРЯ.

2.2.2 Случаи предельных значений параметров.

2.2.3 Нахождение множества блоков фронта

2.2.4 Параллельный алгоритм.

2.2.5 Преимущества построенного метода.

2.3 Алгоритм СРЯ совместно с блочным методом 1ЬИ без решения локальной системы (ВСРЯ).

2.3.1 Замена решения локальной системы предобусловливателем.

2.3.2 Использование единого шаблона для локальной и глобальной матриц.

3 Особенности практической реализации

3.1 Способы хранения матриц.

3.2 Параллелизация алгоритмов.

3.3 Выделение памяти.

3.4 Решение системы линейных уравнений методом РСР

3.5 Решение системы линейных уравнений методом РАСРЫ

3.6 Решение системы линейных уравнений методом РВСРЫ

4 Численные эксперименты

4.1 Гидродинамические модели месторождений.

4.2 Метод РСРЯ.

4.2.1 Поиск оптимального значения £^ос

4.2.2 Масштабируемость.

4.2.3 Тестирование на крупной гидродинамической модели месторождения.

4.2.4 Тестирование совместно с АМС.

4.3 Метод РАСРИ

4.3.1 Сходимость, тестовые гидродинамические модели месторождений.

4.3.2 Комплексное тестирование гидродинамической модели месторождени51.

4.3.3 Случай (1, 1, 1).

4.3.4 Случай (1, 0, 0), влияние начального приближения

4.3.5 Случай (0, 0, 0).

4.3.6 Нахождение оптимального набора параметров

4.3.7 Масштабируемость.

4.3.8 Влияние числа итераций, производимых для расчета фронта.

4.3.9 Влияние скважин на работоспособность алгоритма . 78 4.4 Метод РВСРЯ.

Выводы

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

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

Практически все современные программные пакеты, предназначенные для моделирования описанных процессов, так или иначе сводят исходную задачу к многократному решению системы линейных алгебраических уравнений с сильно разреженной матрицей. При этом возникающие системы имеют очень большую размерность (миллионы неизвестных элементов), а наличие пеоднородностей среды, неструктурированной сетки, искусственных источников и стоков (нагнетательные и добывающие скважины соответственно), гравитационных и капиллярных сил, разломов и трещин еще сильнее усложняет структуру получающихся матриц и ухудшает их алгебраические свойства. Кроме того, постоянное истощение «простых» с точки зрения добычи, а, следовательно, и моделирования, месторождений приводит к необходимости усложнения применяемой математической модели и увеличения степени ее детализации. Следствием этого является то, что большая часть времени и ресурсов ЭВМ часто уходит на решение систем линейных алгебраических уравнений.

Существует множество способов решения систем линейных алгебраических уравнений с разреженной матрицей и их предобусловливания, рассмотренных такими авторами как Г.И. Марчук, Н.С. Бахвалов, Г.М. Кобельков, С.К. Годунов, Р.П. Федоренко, И. Саад и другими. Однако лишь очень немногие из них учитывают специфику данной конкретной задачи (например методы, описанные в работах Д. Уолиса, К. Азиза, К. Штубена). Между тем, именно специфические алгоритмы, использующие различные особенности исходной системы дифференциальных уравнений, практически всегда оказываются наиболее эффективными.

Одним из таких алгоритмов является метод CPR (Constrained Pressure Residual), основанный на идее использования при построении предобусловливателя не только исходной системы линейных алгебраических уравнений, но и матрицы, фактически отвечающей за неявную только по давлению аппроксимацию и потому имеющей в несколько раз меньшую размерность (в 3 раза для трехкомпонентной смеси). Существует множество как академических, так и прикладных публикаций, в которых утверждается, что этот метод, при его использовании совместно с алгебраическим многосеточным алгоритмом, является более эффективным по сравнению с предобусловливателями «общего назначения» (например, неполное LU-разложение или ILU).

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

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

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

Цели работы и задачи исследования

1. Анализ алгоритма СР11 применительно к решению систем линейных алгебраических уравнений, возникающих при моделировании процессов фильтрации вязкой сжимаемой жидкости в пористой среде с учетом факторов, характерных для реальных задача и оказывающих сильное влияние на получающуюся в результате матрицу и, как следствие, на сходимость метода.

2. Разработка на его базе новых, более эффективных с точки зрения скорости работы и универсальности применения алгоритмов.

3. Разработка параллельных программных модулей, реализующих методы.

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

5. Определение оптимальных параметров алгоритмов и степени влияния того или иного свойства пласта и гидродинамической модели на сходимость.

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

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

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

Научная новизна

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

2. Предложены новые, в том числе более эффективные, версии алгоритма, учитывающие, в частности, архитектуру современных

ЭВМ.

3. Разработаны программные модули, реализующие параллельные версии всех рассмотренных методов, которые показали свою применимость при решении реальных задач.

4. Определены преимущества и недостатки алгоритмов, подтвержденные экспериментальным путем при проведении численного моделирования на параллельных ЭВМ с использованием тестовых и реальных гидродинамических моделей газонефтяных пластов. Подтверждена эффективность методов и описан круг задач, для которых использование того или иного способа решения является наиболее оптимальным.

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

Положения, выносимые на защиту

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

2. Разработан подход к программированию, учитывающий особенности предлагаемых методов и архитектуру современных многопроцессорных ЭВМ и позволяющий добиться эффективной работы.

3. Разработаны программные модули, реализующие параллельные версии рассмотренных методов.

4. Выявлены теоретически и подтверждены экспериментально путем численного моделирования на параллельных ЭВМ с использованием тестовых и реальных гидродинамических моделей газопефтяных пластов различные свойства и параметры алгоритмов. Описан круг задач, для которых применение этих методов наиболее оправдано.

Апробация работы. Основные результаты докладывались и обсуждались на:

• научных семинарах каф. вычислительной математики механико-математического факультета МГУ им. М.В. Ломоносова, Москва, 2007-2009 гг.,

• научном семинаре Института вычислительной математики РАН, Москва, 2008 г.,

• научном семинаре кафедры дифференциальных уравнений Московского энергетического института, Москва, 2009 г.,

• конференциях «Ломоносовские чтения», Москва, 2008, 2009 гг.,

• 4-й международной конференции «Математические идеи П.Л. Чебышева и их приложения к современным проблемам естествознания», Обнинск, 2008 г.,

• научном семинаре кафедры математического моделирования МГТУ им. Н.Э.Баумана, Москва, 2009 г.

Публикации. Научные результаты опубликованы в 4 работах, в том числе трех статьях в журналах перечня ведущих рецензируемых научных журналов и изданий ВАК РФ [26, 27, 28]. Журнал [29] вошел в перечень ВАК уже после выхода публикации.

Личный вклад соискателя. Все исследования, изложенные в диссертационной работе и выносимые па защиту, проведены лично соискателем. Во время исследований были использованы визуализаторы, загрузчики и другие функциональные блоки, являющиеся плодом коллективного труда студентов, аспирантов и преподавателей Московского Государственного Университета им. М.В. Ломоносова.

Структура и объем работы. Диссертационная работа состоит из введения, четырех глав, выводов, заключения, приложения и списка литературы. Количество страниц равно 104. Библиография включает 29 наименований.

Выводы

Параллельная версия алгоритма АСРГ1 (и частный случай — СРЯ) демонстрирует свою эффективность на большом числе реальных промышленных задач. Однако, для его применения необходимо знать набор оптимальных параметров, который, как показывают численные эксперименты, меняется при переходе от одной гидродинамической модели к другой. Данный метод может быть полезен в том случае, когда требуется многократно производить одно и то же моделирование, внося небольшие изменения в гидродинамическую модель (добавление скважины, изменение свойств среды, проницаемостей, начальных давлений и так далее). В этом случае можно вначале несколькими запусками подобрать оптимальные значения параметров, а затем использовать их при дальнейших многократных расчетах. Кроме того, уменьшение точек синхронизации и, как следствие, хорошая масштабируемость, позволяют говорить о росте производительности с ростом числа вычислительных ядер.

Что касается параллельного метода ВСРЯ, при правильной реализации с использованием описанных в диссертационной работе оптимизаций он обладает следующими основными преимуществами:

• в отличие от параллельных СРЯ и АСРЯ не имеет регулируемых параметров;

• работает быстрее базового блочного BCGStab-ЬBILU(0) и других рассмотренных алгоритмов;

• стабильность работы такая же, как и у базового блочного алгоритма;

• уменьшение точек синхронизации (в сравнении с базовым) и более равномерное распределение нагрузки между потоками (по сравнению с СРЯ и АСРЯ) положительно сказываются на общей сходимости метода и его масштабируемости;

• размер требуемой памяти незначительно превосходит таковой для блочного метода неполного Ы1 разложения.

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

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

Заключение

Просуммируем основные результаты.

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

2. Разработаны реализации алгоритма СРИ,, адаптированные под архитектуру современных ЭВМ и позволяющие добиться эффективной работы.

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

4. Проведено моделирование на специальных тестовых задачах, показавшее совпадение результатов с известным решением.

5. По результатам проведенного численного моделирования на параллельных ЭВМ с использованием тестовых и реальных гидродинамических моделей газонефтяных пластов определены преимущества и недостатки методов. Описан круг задач, для которых применение того или иного способа решения является наиболее оптимальным.

200 -'-1-'-'

0 20 40 60 80 100

Е1 (ргезэигв %)

40 60 80 100

Е1 (ргевзиге %)

Рис. 4.20: Гидродинамическая модель 5,

20 40 60 80 100 0 20

Е1 (ргеззигв %)

-к

1ос

3 * 10

40 60 80 100

Е1 (ргвввиге %)

Рис. 4.21: Гидродинамическая модель 5, е[ос = 10

Рис. 4.22: Гидродинамическая модель 5, = 5 * 10 3

Рис. 4.35: Гидродинамическая модель 6, £кос = 10 2

1-----1--350

Рис. 4.36: Гидродинамическая модель 6, е\ос = 5*10 3

350

40 60

Е1 (ргеьбиге %)

40 60

Е1 (ргвмигв %)

Рис. 4.37: Гидродинамическая модель 6, = 10 3

1. Wallis J. R. 1.complete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration // SPE 12265, this paper was presented at the Reservoir Simulation Symposium held in San Francisco, CA, November 15-18, 1983

2. Wallis J. R., Kendall R. P., Little Т. E. Constraint Residual Acceleration of Conjugate Residual Methods // SPE 13536, this paper was presented at the 1985 Reservoir Simulation Symposium held in Dallas, Texas, February 10-13, 1985.

3. Stuben К. Algebraic Multigrid (AMG): Experiences and Comparisons // proceedings of the International Multigrid Conference, Copper Mountain, CO, April 6-8, 1983.

4. Christian Wagner Introduction to Algebraic Multigrid // Course Notes of an Algebraic Multigrid Course at the University of Heidelberg in the Wintersemester 1998/99 Version 1.1.

5. Yousef Saad Iterative Methods for Sparse Linear Systems. Second edition with corrections. January 3rd, 2000.

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

7. Klaus Stueben, D-5205 St. Augustin (W.-Germany), John Rüge, Fort Collins (USA), Institute for Computational Studies at CSU, Rolf Hempel, D-5205 St. Augustin (W.-Germany), AMG1R5 function fortran code.

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

9. Stuben K. Algebraic Multigrid (AMG): an introduction with applications 11 GMD Report 70, 1999.

10. Горелов И. Г. Параллельный адаптивный предобусловливатель СРЯ в задаче фильтрации // Программные продукты и системы, номер 3, 2009 г, стр. 107-109.

11. Богачев К.Ю, Горелов И. Г. Применение параллельного предобусловливателя СРЯ к задаче фильтрации вязкой сжимаемой жидкости в пористой среде / / Вычислительные методы и программирование, 2008, том 9, стр. 184-190.