CC BY
214
УПРАВЛЕНИЕ, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И ИНФОРМАТИКА
УДК 519.6:622.32
А. А. Яппарова, С. А. Маякова
РАСПАРАЛЛЕЛИВАНИЕ АЛГОРИТМОВ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ПРОЦЕССОВ ПЕРКОЛЯЦИИ С ВЫТЕСНЕНИЕМ
Рассматриваются методы параллельной реализации алгоритмов численного моделирования процессов перколяции с вытеснением. Исследуется влияние граничных условий на фрактальную размерность перколяционного кластера. В зависимости от типа граничных условий выбираются различные подходы к распараллеливанию расчетных программ: распараллеливание по статистике и по пространству. Проводится сравнение ускорений, полученных для каждого из методов распараллеливания. Двухфазная фильтрация; перколяция с вытеснением; фрактальная размерность кластера; параллельные алгоритмы
ВВЕДЕНИЕ
В различных прикладных задачах нефтяной промышленности производится моделирование течения флюидов (нефти, воды) в пласте. Существует множество моделей, использующих различные подходы к описанию физической сущности процессов, протекающих в нефтяном пласте. Модель перколяции с вытеснением с образованием захваченных областей достаточно хорошо отражает физические свойства процесса фильтрации и позволяет прослеживать движение фронта раздела вытесняющей и вытесняемой жидкостей. Кроме того, возможность образования «захваченных» областей позволяет учитывать несжимаемость нефти и оценивать остаточную нефтенасыщенность.
Рассмотрим вытеснение нефти водой. Если вода нагнетается очень медленно, то капиллярные силы преобладают над вязкими и динамика вытеснения полностью определяется процессами в масштабе отдельных пор. В пределе можно пренебречь всеми скачками давления как в вытесняющей (воде), так и в вытесняемой жидкости (нефти). Однако между этими двумя жидко -стями сохраняется разность давлений (капиллярное давление), равное
Р Р _ 2о ^0 (1)
Р вытесняющее Р вытесняемое г ’ ' '
здесь о - межфазное натяжение на границе жидкостей, 0 - краевой угол между поверхностью раздела и стенкой поры, а г - радиус кривизны поры в месте контакта жидкостей.
Капиллярные силы наиболее велики в самых узких местах системы пор. Таким образом, если горловины уже пор, то поверхность раздела во-
Контактная информация: (347) 273-36-22 Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии-2010»
ды и нефти быстро движется сквозь горловины и тормозится, попадая в крупные поры. В согласии с подобными простыми теоретическими моделями и экспериментальными результатами это движение можно представить в виде ряда дискретных скачков, так что на каждом временном шаге вода вытесняет нефть из самой малой из имеющихся пор, занятых нефтью.
Вытесняющая жидкость захватывает области вытесняемой жидкости. По мере продвижения вторгающейся жидкости она способна полностью охватить области, заполненные вытесняемой средой, т. е. полностью отрезать кластеры вытесняемой жидкости от узлов на входе в решетку. Это одна из причин возникновения проблемы «остаточной нефти» - жгучей экономической проблемы нефтепромыслов. Поскольку нефть несжимаема, можно сформулировать правило: вода не может вытеснить нефть из захваченных областей.
Для описания течения двух несмешиваю-щихся жидкостей в пористой среде (в нашем случае это процесс вытеснения нефти водой) Уилкинсон и Виллемсен в начале 1980-х годов предложили модель перколяции с вытеснением (invasion percolation) [1-5]. Уилкинсон и Вил-лемсен предложили моделировать этот процесс в идеализированной среде, в которой пористую среду можно рассматривать как регулярную решетку с узлами и связями, играющими роль пор и связывающих их каналов (горловин). Физическая неоднородность привносится в эту среду путем приписывания узлам и связям чисел, которые характеризуют размеры соответствующих пор и горловин. Моделирование процесса на выбранной реализации решетки состоит, таким образом, в расчете движения поверхности раздела воды и нефти по мере того, как она продвигается через самые малые из имеющихся пор, заполняя их вытесняющей жидкостью. В модели перколяции с вытеснением
и возможностью образования захваченных областей (trapping invasion percolation) учитывается практическая несжимаемость нефти, а именно: области, полностью окруженные водой, считаются «захваченными», и нефть не может быть извлечена из таких областей.
Перколяция с вытеснением - динамический перколяционный процесс, требующий пошагового расчета роста кластера. Кластер вытесняющей жидкости, полученный в результате моделирования процесса, имеет фрактальную структуру. Для оценки различных фрактальных свойств перколяционного кластера (например, фрактальной размерности) необходимо абстрагироваться от конкретной реализации решетки (вместо конкретных размеров пор и горловин, узлам и связям решетки присваиваются случайные числа) и провести большое количество испытаний на различных реализациях решетки с последующим усреднением. Среди методов численного решения задачи о перколяции особое место занимают методы Монте-Карло, получившие широкое распространение в связи с увеличивающимися темпами развития вычислительных мощностей.
1. МОДЕЛЬ ПЕРКОЛЯЦИИ С ВЫТЕСНЕНИЕМ
Рассматривается процесс перколяции с вытеснением на решетке размером 2L х 3L. Вытесняющая жидкость впрыскивается на левой границе решетки, процесс прослеживается до тех пор, пока не наступает протекание, т. е. пока вытесняющая жидкость не достигнет правой границы (рис. 1).
Рис. 1. Модель эксперимента
Исследуется влияние непроницаемой нижней границы на поведение растущего кластера. Верхняя и нижняя граница оказывают аналогичное влияние на кластер, поэтому, чтобы исключить влияние верхней границы, вытесняющая жидкость впрыскивается только в нижние Ь ячеек левой границы. Изучаются свойства левого нижнего квадрата размера Ь х Ь, в котором
можно считать, что процесс установился и асимптотические фрактальные особенности перколяционного кластера ярко выражены.
2. ПОСЛЕДОВАТЕЛЬНЫЙ АЛГОРИТМ
Обратимся к подробному описанию моделирования процесса перколяции с вытеснением. Рассматривается двумерная регулярная прямоугольная решетка. Каждой ячейке решетки ставится в соответствие число г из интервала (0, 1). Формируется двумерный массив, значение элемента массива равно нулю, если соответствующая ему ячейка заполнена нефтью. Если в ячейке вода, то значение соответствующего элемента массива равно положительному числу, равному номеру шага, на котором произошло вытеснение нефти водой из соответствующей ячейки. Выберем место впрыскивания вытесняющей жидкости (воды) - источник - и место вытекания вытесняемой среды (нефти) - сток. Элементам массива, соответствующим местам впрыскивания, присвоим значения, равные единице.
Алгоритм
Найдем узлы роста как узлы, занятые вытесняемой жидкостью и соседствующие узлам с вытесняющей средой.
Пропустим вытесняющую среду в тот узел роста, в котором случайное число г принимает наименьшее значение. Присвоим этому узлу значение, равное номеру шага.
При моделировании перколяции с захватом на каждом шаге необходимо проверять наличие захваченных областей. Области, полностью окруженные вытесняющей жидкостью, помечаются как «захваченные» и исключаются из списка узлов роста, для обнаружения таких областей используется алгоритм Хошена-Копельмана.
Возвращаемся ко второму шагу алгоритма и повторяем процесс до тех пор, пока вытесняющая жидкость не достигнет стока.
Поясним подробнее пункт 3 алгоритма. Для того чтобы определить «захваченные» области, необходимо найти на решетке все кластеры нефти и исключить кластеры, соединенные с открытой границей. Все оставшиеся кластеры нефти будут со всех сторон окружены водой, и нефть из них вытеснить будет уже никак нельзя. Поиск кластеров нефти проводится по алгоритму Хошена-Копельмана [6].
Важно отметить, что принадлежность узла какому-либо кластеру - глобальное свойство и может быть определена только после просмотра всей решетки, а главное достоинство
алгоритма Хошена-Копельмана в том, что он позволяет идентифицировать все кластеры за один проход по решетке. Идея алгоритма состоит в том, что узлам, принадлежащим к одному кластеру, присваиваются метки с одинаковым номером.
3. АЛГОРИТМ ХОШЕНА-КОПЕЛЬМАНА МАРКИРОВКИ КЛАСТЕРОВ
Последовательно обходим все ячейки решетки по строкам, начиная с первой. Будем называть ячейку свободной, если в ней находится вытесняемая жидкость, и занятой, если в ней находится вытесняющая жидкость. Производится поиск и маркировка свободных ячеек. Необходимо ввести дополнительный массив «правильных» меток, к-й элемент этого массива будет равен к, если к-я метка - «правильная», и равен г (г < к), если все ячейки, помеченные номером к, на самом деле принадлежат кластеру с меткой г. Начинаем присваивать ячейкам кластерные метки, начиная с единицы. Для каждой свободной ячейки существует четыре возможных варианта:
Соседняя ячейка слева занята, соседняя сверху занята (в качестве рабочей гипотезы принимаем, что данная ячейка входит в новый кластер, присваиваем ячейке текущий номер кластерной метки, увеличиваем номер кластерной метки на единицу; конечно, может оказаться, что новый, как нам кажется, кластер является просто ответвлением какого-то другого кластера, мы узнаем об этом, только просмотрев всю решетку).
Соседняя ячейка слева свободна, соседняя сверху занята (текущая ячейка и ее сосед слева принадлежат к одному кластеру; присваиваем ячейке метку соседа слева).
Соседняя ячейка слева занята, соседняя сверху свободна (текущая ячейка и ее сосед сверху принадлежат к одному кластеру, однако у соседа сверху кластерная метка может быть «неправильной», так как в результате проверок могло выясниться, что кластер, к которому принадлежит эта ячейка, слился с другим кластером; присваиваем текущей ячейке метку соседа сверху, либо «правильную» метку соседа сверху).
Соседняя ячейка слева свободна, соседняя справа свободна (если метки соседей справа и слева не совпадают, то считаем большую из них «неправильной» и заносим в соответствующий элемент массива «правильных» меток номер меньшей из меток соседей; присваиваем
текущей ячейке метку, наименьшую из «правильных меток» соседних справа и слева ячеек).
4. ФРАКТАЛЬНАЯ РАЗМЕРНОСТЬ КЛАСТЕРА
На рис. 2 представлена одна из реализаций эксперимента при Ь = 40. Из рисунка видно, что перколяционный кластер заполняет некоторую площадь, но имеет сложную структуру: сильно изрезанный внешний периметр, множество «дырок» во внутренней части.
Рис. 2. Типичная реализация эксперимента, Ь = 40
Перколяционный кластер обладает фрактальной размерностью. Число занятых узлов решетки в момент наступления протекания растет с размером решетки по закону
Ы(Ь) = а ■ , (2)
где й - фрактальная размерность кластера при перколяции с вытеснением с образованием захваченных областей. В различных работах, посвященных решению задач о перколяции, приводится следующая оценка значения фрактальной размерности кластера - й = 1,82 ± 0,01 [2, 3].
Изучая фрактальную структуру кластера, мы ищем легкий путь описать ее основные свойства. Наиболее простой способ - выбрать подмножество меньшей размерности, которое полностью характеризует всю структуру. Для фрактала размерности В, вложенного в двумерную область, это точки пересечения структуры и прямой. Подмножество пересечения также фрактально, и его размерность й = В - 1.
В данной работе изучалась фрактальная размерность набора точек, полученных пересечением кластера линией, параллельной нижней границе. Прослеживается изменение фрактальной размерности кластера от границы к средней области.
5. АНАЛИЗ РЕЗУЛЬТАТОВ
В ходе численного эксперимента первой задачей было оценить массовую размерность кластера. Для этого велся подсчет ячеек, занятых вытесняющей жидкостью. Причем статистика собиралась только по левому нижнему квадрату размера Ь х Ь, так как можно считать, что в этой области процесс установившийся.
Однако основной целью работы была проверка гипотезы об уменьшении фрактальной размерности кластера в направлении от средней области решетки к непроницаемой границе [7, 8]. Интуитивно понятно, что вероятность роста кластера около границы меньше, так как у граничных ячеек на одного соседа меньше, чем у центральных, то есть около границы меньше узлов роста, а раз кластер, вероятнее всего, будет расти в центре решетки, то размерность его около границы должна быть меньше.
Для повышения достоверности экспериментальные данные осреднялись по статистике. Была проведена серия из 1000 численных экспериментов для значений Ь = 32, 64, 128.
При расчетах было проведено сравнение результатов, полученных с применением стандартного генератора случайных чисел и генератора МТ19937-64 на основе алгоритма твистера Мерсенна. Оценки фрактальных размерностей, полученные с помощью генератора МТ19937-64, лучше согласуются с теоретическими.
Теперь обратимся к анализу численных результатов. Прологарифмировав (2), получим
1пЩ (Ь) = й1пЬ + 1па, (3)
т. е. фрактальная размерность кластера й - это угловой коэффициент прямой, выражающей зависимость Щ(Ь) в двойных логарифмических координатах (см. рис. 3).
На рис. 3 приведено уравнение прямой, линейно аппроксимирующей искомую зависимость. Была получена оценка фрактальной размерности методом подсчета клеток й~ 1,8174, что вполне согласуется с экспериментально полученным значением й =1,82 из различной литературы, посвященной перколяции с вытеснением с образованием захваченных областей.
Одной из важных задач, продиктованных целью работы, было проследить влияние границы на размерность кластера, а именно рост фрактальной размерности пересечения при движении от нижней границы к внутренней области кластера. Это явление было прослежено в ходе эксперимента, результаты наглядно представлены на рис. 4, здесь г изменяется от единицы до Ь, таким образом, г / Ь - нормализованное расстояние от границы. Исследовалось
поведение величины 1п Щ(Ь) / 1п Ь, близкой к значению размерности с точностью до 1п а / 1п Ь ^ 0 при Ь ^да.
1п N(1) у= 1,8174*-0,3505
Рис. 3. Зависимость числа ячеек перколяционного кластера от размера решетки в двойных логарифмических координатах. Массовая размерность кластера оценивается как угловой коэффициент наклона прямой
Рис. 4. Изменение фрактальной размерности пересечения в зависимости от расстояния от границы
Важным доказательством адекватности полученных результатов является то, что верхняя граница изменения размерности пересечения согласуется с полученным значением массовой размерности, т.е. фактически экспериментальное значение массовой размерности не зависит от способа измерения.
6. ПОДХОДЫ К РАСПАРАЛЛЕЛИВАНИЮ
При моделировании процессов методом Монте-Карло необходимо проведение большого числа статистических испытаний. Естественно, с увеличением числа испытаний, время работы программы линейно растет. Совершенно очевидной была идея произвести распараллеливание нашей программы по испытаниям. Это было
сделано средствами ОрепМР с минимальными изменениями в коде программы.
Распараллеливание по статистике, как и следовало ожидать, дает очень хорошие результаты по ускорению и эффективности (см. табл. 1). Очевидно, что распараллеливание по статистике дает линейный рост ускорения, что было проверено на четырехядерной системе и должно выполняться и для систем с большим количеством ядер.
Т аблица 1
Число ядер Время работы программы, с Ускорение Эффективность
1 821,78 - -
2 414,59 1,98 0,99
3 278,99 2,95 0,98
4 211,08 3,89 0,97
Хотя для большей достоверности без проведения большого числа испытаний не обойтись, важно отметить, что оценки фрактальных размерностей тем точнее, чем больше размер решетки, так как размерность Хаусдорфа вычисляется в нашем случае по формуле d = N (L)
= lim------. Но с увеличением размера решет-
l
ки, время работы программы растет нелинейно (см. рис. 5), и, начиная с некоторых значений L, мы столкнулись с тем, что время расчета одного испытания на одном процессоре слишком велико.
Все приведенные ниже результаты были получены на вычислительной платформе со следующими характеристиками: 2 x Intel Xeon 2.00 GHz (4 ядра), 4Mb Cache, 4Gb RAM.
Рис. 5. Динамика роста времени расчетов одного испытания на одноядерном процессоре с увеличением размера решетки
В связи с этим мы стали рассматривать различные варианты пространственного распараллеливания по решетке каждого отдельного ис-
пытания. После тщательного анализа работы программы было выяснено, что наибольшее время счета занимает функция поиска «захваченных» областей. Таким образом, проблема распараллеливания отдельного эксперимента была решена использованием на каждом шаге роста кластера не последовательной, а параллельной реализации алгоритма Хошена-Копельмана [9].
7. ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ АЛГОРИТМА ХОШЕНА-КОПЕЛЬМАНА
Решетка делится на полосы в направлении роста кластера, и каждый процессор исследует только свою полосу. Маркировка кластеров на каждом процессоре выполняется локально и только в конце необходим обмен граничными данными между процессорами.
Кластеры, полностью расположенные внутри одной полосы и не соединяющиеся с границами, естественно называть локальными. Кластеры, расположенные на нескольких полосах и хранящиеся соответственно на нескольких процессорах, будем называть глобальными.
В алгоритме Хошена-Копельмана каждому кластеру присваивается метка с номером. Различают «правильные» и «неправильные» метки. «Неправильные» метки появляются в процессе подсчета, когда кажется, что обнаружен новый кластер, а оказывается, что это часть уже известного кластера. Для «неправильных» меток известен номер соответствующей «правильной» метки.
Для параллельной реализации необходимо ввести еще один тип меток - «глобальные». В этих структурах данных должна храниться информация о том, какие локальные кластеры связаны между собой, то есть образуют глобальный кластер.
В табл. 2 приведены результаты пространственного распараллеливания.
Таблица 2
Число ядер Время работы программы, с Ускорение Эффективность
1 5412,73 - -
2 3689,32 1,47 0,74
3 2700,66 2,00 0,67
4 2183,87 2,49 0,62
Ускорение и эффективность конечно меньше, чем при распараллеливании по статистике, потому что при пространственном распаралле-
ливании необходимо наличие последовательных участков кода. Эффективность падает линейно в связи с линейным ростом количества областей сшивки кластера по закону (2p — 3), где р - число ядер. Получив линейную зависимость для четырехядерной системы, мы можем сделать вывод о линейном характере роста ускорения и для большего количества ядер.
ВЫВОДЫ
Оба рассмотренных подхода к распараллеливанию моделирования процесса перколяции с вытеснением доказали на практике свою эффективность и полезность. Анализ ускорения и эффективности показал, что распараллеливание по числу статистических испытаний выгодно осуществлять, когда NMC / L > 10, где NMC - число экспериментов. Начиная с L = 256, когда условие NMC / L > 10 нарушается, распараллеливание только по статистике уже невыгодно, так как расчет одного эксперимента занимает продолжительное время. При моделировании процесса на больших размерах решеток необходимо пространственное распараллеливание. Надо отметить, что возможно успешное комбинирование этих двух подходов при больших размерах решеток и большом количестве статистических испытаний.
СПИСОК ЛИТЕРАТУРЫ
1. Percolation and capillary fluid displacement / Koplik J. [and etc] // Schlumberger-Doll Research preprint, 1983. P. 1-15.
2. Sahimi M. Applications of percolation theory. Taylor&Francis, 1993. 258 p.
3. Stauffer D., Aharony A. Introduction to percolation theory. Taylor&Francis, 2003. 181 p.
4. Федер Е. Фракталы. М.: Мир, 1991. 254 с.
5. Шредер М. Фракталы, хаос, степенные законы. Миниатюры из бесконечного рая. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001. 528 с.
6. Тарасевич Ю. Ю. Перколяция: теория, приложения, алгоритмы: Учеб. пособие. М.: Едиториал УРСС, 2002. 112 с.
7. Surface effects in invasion percolation / Cafie-ro R. [and etc] // Phys. Rev. E. 1997. P. 1-4.
8. Theory of boundary effects in invasion percolation / Gabrielli A. [end etc] // arxiv.org. 1998. P. 1-11.
9. Tiggemann D. Simulation of percolation on mas-sively-parallel computers // arxiv.org. 2008. P. 1-9.
ОБ АВТОРАХ
Яппарова Алина Анваровна,
студ. каф. высокопроизв. вы-числ. технол. и систем. Иссл. в обл. матем. моделирования.
Маякова Светлана Алексеевна, доц. той же каф. Дипл. ма-тем., системн. программист
(УГАТУ, 2004). Канд. физ.-мат. наук по матем. моделированию, числ. методам и комплексам программ (УГАТУ, 2008). Иссл. в обл. матем. моделирования.
