CC BY
14
Вычислительные технологии
Том 22, № 3, 2017
Высокопроизводительные гетерогенные вычисления CPU — GPU в задаче электрического каротажа нефтегазовых скважин
В.Н. Глинских1'2'*, А. Р. Дудаев1'2, О. В. Нечаев1
1Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, Россия
2Новосибирский национальный исследовательский государственный университет, Россия *Контактный e-mail: glinskikhvn@ipgg.sbras.ru
Разработан и программно реализован алгоритм моделирования данных электрического каротажа нефтегазовых скважин с использованием высокопроизводительных гетерогенных вычислений CPU-GPU. Численное решение прямой двумерной задачи основано на методах конечных элементов и разложении Холецкого для решения системы линейных алгебраических уравнений (СЛАУ). Программные реализации алгоритма созданы с использованием технологии NVIDIA CUDA и вычислительных библиотек, позволяющих выполнять разложение СЛАУ и находить ее решение на CPU и GPU. Проведен анализ времени вычислений в зависимости от размера матрицы и количества ее ненулевых элементов. Получены оценки быстродействия вычислений на CPU и GPU, а также высокопроизводительных гетерогенных CPU-GPU. С использованием разработанного алгоритма проведено моделирование данных электрического каротажа в реалистичных геологических средах.
Ключевые слова: графический процессор, параллельный алгоритм, гетерогенные вычисления CPU-GPU, конечно-элементный метод, прямая двумерная задача, электрический каротаж.
Введение
В последнее десятилетие значительно усложняется и расширяется круг задач сква-жинной электрометрии. Изучаются типы залежей нефти и газа, имеющих сложное геологическое строение. Это требует совершенствования каротажных приборов и их программно-алгоритмического обеспечения, необходимого для обработки и интерпретации скважинных данных. Одной из каротажных технологий, широко используемых при изучении нефтяных и газовых скважин, является аппаратурный комплекс СКЛ [1]. С его применением каротаж выполняется за одну спускоподъемную операцию, что сокращает время задалживания скважины. В едином комплексе реализуются основные геофизические методы исследований в скважинах, среди которых особое место занимает электрический каротаж.
Внедрение новой аппаратуры приводит к необходимости разработки современных программ и алгоритмов с целью повышения эффективности применения каротажных методов и точности интерпретации данных. Особо важно увеличение точности оценки
© ИВТ СО РАН. 2017
нефтегазонасыщения пластов, которая выполняется на основе данных электрического каротажа по значениям удельного электрического сопротивления (УЭС) горных пород. Для этого требуется использование многомерного численного моделирования в рамках реалистичных моделей геологических сред на основе решений задач электродинамики в полных математических постановках.
Для решения вычислительных задач электрического каротажа используются различные методы. Широкое распространение имеют численно-аналитические и приближенные методы, поскольку позволяют создать на своей основе быстрые алгоритмы обработки данных [2-4]. Однако при моделировании электрических полей в пространственно неоднородных средах основными являются конечно-разностный и конечно-элементный методы [5,6]. Как известно, численные решения задач электродинамики с использованием сеточных методов приводят к решению систем линейных алгебраических уравнений (СЛАУ) с разреженными матрицами большой размерности, для чего эффективно применяются прямые и итерационные методы. Но использование этих решений для оперативной интерпретации данных в масштабе реального времени оказывается слабоэффективным в силу их низкой производительности и высокой ресурсоемкости.
Одним из способов сокращения времени расчетов является применение многопроцессорных компьютерных систем и вычислительных кластеров. Их использование для оперативного решения интерпретационных задач каротажа [7], в том числе непосредственно на скважине, достаточно затруднительно [8]. Однако при наличии соответствующих каналов связи применение вычислительных кластеров возможно и в удаленных районах. В этом случае решение задачи переносится на вычислительные системы более высокого класса по сравнению с рассмотренными в статье. По этой теме существует большое число публикаций, а полный обзор всего массива работ по вычислениям на кластерах мог бы составить отдельную книгу. Решениям задач прикладной геофизики, в том числе промысловой, с использованием многопроцессорных вычислительных систем, включая системы Сибирского компьютерного центра СО РАН, посвящены работы [6,9-12].
Другой способ связан с вовлечением графических процессоров NVIDIA и сопроцессоров Intel Xeon Phi для ускорения вычислений. Он широко применяется при решении современных задач в разных научных областях, в том числе в каротаже [13-17]. Основная идея представленной работы заключается в изучении возможностей использования мобильных многопроцессорных устройств (графических процессоров) применительно к важным практическим задачам каротажа, решение которых необходимо в масштабе реального времени, в первую очередь в процессе бурения нефтяных и газовых скважин.
Повышение эффективности применения численных решений задач электродинамики обусловлено необходимостью создания высокоэффективных параллельных алгоритмов. Последние в большей степени связаны с разработкой или применением быстрых вычислительных процедур и библиотек для решения СЛАУ, основанных на различных методах. Известно решение рассматриваемой в работе задачи методом конечных разностей [16]. Задача сводится к решению СЛАУ с положительно определенной матрицей и после симметризации эффективно решается методом сопряженных градиентов с предобусловливанием на основе алгоритма Шульца — Хотеллинга. Реализация этого алгоритма на GPU позволила сократить время вычислений более чем на один порядок по сравнению с последовательным алгоритмом. Другие результаты успешного применения GPU для решения задач электрокаротажного зондирования в скважинах приведены, например, в [18,19].
В настоящей работе для нахождения решения СЛАУ используется прямой метод на основе разложения Холецкого и последующего решения двух вспомогательных СЛАУ с треугольными матрицами. Традиционно при разработке таких высокопроизводительных параллельных прямых алгоритмов большое внимание уделяется эффективности разложения матрицы. Для этого разрабатываются как сами вычислительные алгоритмы, так и соответствующие структуры данных [20]. Это связано с тем, что решение вспомогательных СЛАУ занимает гораздо меньше вычислительных ресурсов по сравнению с разложением основной матрицы. Разработанные параллельные алгоритмы разложения матриц реализуются также и на графических процессорах (GPU) [21,22]. Но в отличие от центральных процессоров (CPU), GPU позволяют эффективно реализовать довольно узкий класс алгоритмов, что приводит к необходимости создания методов, использующих одновременно CPU и GPU [23]. В работах [20,24] рассматриваются подобные гетерогенные алгоритмы разложения матриц.
В случае, когда необходимо решать много СЛАУ с одной и той же матрицей, но с разными правыми частями, вычислительные затраты на решение двух СЛАУ с треугольными матрицами значительно превышают затраты на разложение матрицы основной СЛАУ, что может приводить к низкой эффективности решений таких СЛАУ при использовании только GPU. Отсюда возникает необходимость разработки гетерогенного метода решения подобных СЛАУ, использующего вычисления CPU-GPU.
В конечном счете эти подходы делают возможным применение решений ресурсоемких вычислительных задач для практических целей обработки данных электрического каротажа, обеспечивая повышение точности определения параметров нефтегазового пласта. Представленная работа посвящена разработке алгоритма и его программной реализации решения прямой двумерной задачи электрического каротажа нефтегазовых скважин на основе метода конечных элементов и высокопроизводительных гетерогенных вычислений CPU - GPU.
1. Решение прямой двумерной задачи БКЗ
При изучении геологических разрезов большую роль играет определение УЭС горных пород вокруг скважины для оценки их нефтегазонасыщения. Один из самых известных способов решения этой задачи состоит в регистрации электрических потенциалов, наводимых источником постоянного тока в скважине. Среди методов электрического каротажа особое место занимает боковое каротажное зондирование (БКЗ) [25]. В этом случае скважинный прибор состоит из нескольких соосно размещенных зондов: одного питающего (токового) электрода А и двух сближенных измерительных М- и Ж-электродов. Использование зондов разной длины обеспечивает изучение среды в радиальном направлении. Такое зондирование необходимо для получения распределения УЭС от стенки скважины в глубь пласта с последующим определением характера флюидонасыщения горных пород.
В методе БКЗ кажущееся удельное электрическое сопротивление вычисляется по формуле
Рк = к-—-, (1)
где к — геометрический коэффициент зонда, определяемый расстояниями между электродами; ^м — потенциал, измеренный на электроде М; — потенциал, измеренный на электроде N; — сила тока, протекающего через токовый электрод А.
А
М N
о
■о-о
2.0 м
0.5 м
Рис. 1. Схема расположения зондов БКЗ
Результатом измерений БКЗ являются каротажные диаграммы, отражающие зависимость значений кажущегося УЭС от глубины скважины. В БКЗ используются разноглубинные трехэлектродные зонды, которые имеют обозначение A2.0M0.5N. Здесь расстояние между электродами А и М составляет 2.0м, а между М и N — 0.5м (рис. 1). В работе приводится анализ диаграмм, полученных при использовании зондов длиной от 0.4 до 8.0 м.
Основные идеи и подходы физических измерений, выполняемых методом бокового каротажного зондирования, заложены еще в прошлом веке, но их дальнейшее развитие во многом зависит от существующей математической базы. Учет эффектов, главным образом связанных со сложным строением геологической среды, требует вовлечения многомерного математического аппарата.
Рассмотрим решение прямой двумерной задачи БКЗ. Двумерная осесимметричная геоэлектрическая модель среды описывает геологический разрез, вскрытый вертикальной скважиной (рис. 2). Модель включает пачку пластов с плоскопараллельными горизонтальными границами, пересеченными вертикальной цилиндрической скважиной. В прискважинной области каждого пласта может быть зона проникновения, сформированная в процессе вытеснения пластового флюида фильтратом бурового раствора. Она отделена от скважины и пласта коаксиально-цилиндрической границей. Каждая из областей геоэлектрической модели среды характеризуется своими значениями УЭС.
Традиционно под решением прямой задачи БКЗ или моделированием каротажных измерений понимают определение кажущегося УЭС по заданной функции пространственного распределения УЭС среды, описываемого фиксированным вектором модельных параметров. Как следует из уравнения (1), для нахождения кажущегося УЭС необходимо знать электрические потенциалы на измерительных электродах М и N.
Рис. 2. Геоэлектрическая модель среды с осевой симметрией, описывающая флюидонасыщен-ный пласт-коллектор с проникновением фильтрата бурового раствора
Для нахождения электрического потенциала в области моделирования решается следующая краевая задача:
—а grad = 0, (2)
а
а
др д п
др д п
0,
(3)
г1
0,
3,
(5)
г2
где — потенциал напряженности электрического поля Е = — grad р; а = 1/р — удельная электрическая проводимость; Г0 — внешняя граница области, на которой электрический потенциал считается близким к нулю; Г1 — диэлектрическая поверхность зонда; Г2 — поверхность токового электрода; ] — плотность тока.
Поскольку рассматриваемая осесимметричная модель имеет цилиндрические границы, удобно от декартовой системы координат перейти к цилиндрической, тогда дифференциальные операторы примут следующий вид:
grad и =
ди 1 ди ди дг' г 8ф' дг
)
1 дгУг 1 дУф dlv V =---1---— Н--.
^^ ¥ * г\ 1 г\ / 1 О
г дг г дф дг
Далее полагаем, что область моделирования обладает осевой симметрией, а оси скважины и зонда совпадают, тогда краевая задача (2)-(5) примет вид
1 д / д /
г дг \ дг ) дг \ дг)
(Ун
^|Г0 = 0,
др дг
(6) (7)
0,
г1
др дг
= 3.
г2
Введем функциональные пространства:
Н 1(П) = {и е Ь2(П) | grad и е Ь2(П)} ,
#0(П) = {и е Н 1(П) | м|го = 0} . Для элементов этих пространств введем следующее скалярное произведение:
(и,ь)
ии ¿П.
Для решения краевой задачи (6)-(9) используем метод конечных элементов [26]. Выражение (6) умножим скалярно на некоторую функцию V е Н^ (П) и, воспользовавшись формулой интегрирования по частям, а также тем фактом, что функция V равна нулю
на части границы расчетной области, в результате получим следующую вариационную постановку для краевой задачи (6)-(9).
Найти tp Е Щ (П) такое, что V v Е Щ (П) выполняется соотношение
1° (II+Ш)м (1°)
п г2
Приближенное решение в вариационной постановке (1°) будем искать в виде разложения по множеству базисных функций фг, образующему конечномерное подпространство пространства Н^П). В качестве базиса используем билинейные базисные функции [26]. В результате приближенное решение полностью определяется состоящим из весов разложения вектором x, найти который можно, решив следующую СЛАУ:
Ax = b.
Элементы матрицы и вектора правой части можно вычислить как
mI = ia (дф1дф1 + дф1дф1 \ rdrdz [A]« = J а\дг дг + dz dz )rdrdZ' п
[Ь]г = J JAФг dT2. г2
Для решения результирующей СЛАУ используем разложение Холецкого
A = LL\ (11)
где L и LL — верхняя и нижняя треугольные матрицы.
Необходимо отметить, что при движении зонда вдоль скважины во всей краевой задаче (6)-(9) будет меняться лишь расположение границ Г:и Г2, а это, в свою очередь, ведет к изменению правой части уравнения и неизменности элементов матрицы A. Следовательно, для моделирования измерений в нескольких точках скважины разложение матрицы на произведение двух треугольных матриц можно выполнить один раз, а затем использовать его для нахождения нескольких решений СЛАУ.
Таким образом, решение прямой задачи БКЗ позволяет по заданному вектору параметров геоэлектрической модели вычислить значения кажущегося УЭС на заданной глубине вдоль скважины. Результатом решения соответствующего количества прямых задач являются синтетические данные, представляемые в виде диаграммы.
2. Оценка быстродействия вычислений на CPU и GPU в задаче БКЗ
На основе специализированных вычислений на GPU персональных компьютеров разработаны параллельные алгоритмы моделирования данных БКЗ в нефтегазовых скважинах. С их помощью получены оценки быстродействия вычислений на CPU (Intel Core i7-382°) и GPU (NVIDIA GeForce Titan), включая гетерогенные CPU-GPU. Используемые CPU и GPU имеют следующие характеристики: число физических ядер 4,
виртуальных ядер 8, частота процессора 3.6 ГГц; число CUDA-ядер 2688, частота процессора 876 МГц.
Разработка параллельных алгоритмов выполнена с использованием технологии NVIDIA CUDA [27] и вычислительных библиотек cuSOLVER и CHOLMOD. Последние имеют расширенный интерфейс программирования приложений API (Application Programming Interface). В библиотеке cuSOLVER количество вычислительных потоков на CPU и GPU может устанавливаться автоматически в зависимости от вычислительных операций и требуемых для этого ресурсов. При работе с библиотекой CHOLMOD на CPU выбирается стандартное количество вычислительных потоков (10), а на GPU — максимально возможное (2688).
Далее проводится анализ времени вычислений при решении одной прямой задачи БКЗ в заданной точке по профилю скважины для одного зонда. При этом с использованием (11) выполняется факторизация — разложение матрицы СЛАУ на треугольные матрицы L и L*, которое требуется для нахождения одного решения СЛАУ. Поскольку основные вычислительные затраты связаны с факторизацией и составляют более 80 % от общего времени исполнения программы, эти вычисления целесообразно выполнять параллельно на GPU.
Разработан параллельный алгоритм моделирования данных БКЗ на основе вычислительной библиотеки cuSOLVER, интегрированной в состав CUDA Toolkit 7 и выше [28]. С ее применением осуществляются как факторизация, так и нахождение одного или несколько решений СЛАУ непосредственно на GPU. При этом версия cuSOLVER, позволяющая выполнять их раздельный вызов, представлена в CUDA Toolkit 7.5. Использование разработанного алгоритма показало, что в случае решения одной прямой задачи БКЗ время вычислений при факторизации на GPU сокращается более чем в два раза.
С помощью NVIDIA Visual Profiler выполнен анализ соотношения временных затрат при вычислениях на GPU (рис. 3). Так, время передачи данных с Host на Device
Рис. 3. Соотношение временных затрат при вычислениях на GPU: 1 — передача данных с Host на Device и обратно; 2, 3 — анализ треугольных матриц L и L1; 4 — сортировка элементов; 5 — факторизация; 6, 7 — использование треугольных матриц L и L1 при нахождении решения СЛАУ
и обратно составляет около 1 % от всех операций на GPU. Анализ треугольных матриц L и Lt и сортировка элементов составляют от 3.5-4.5 % от общего времени. Наибольшее время занимает факторизация — около 70 %. Нахождение решения СЛАУ с использованием двух треугольных матриц требует около 18 % в рассматриваемом вычислительном эксперименте при решении одной прямой задачи БКЗ.
Как видно из проведенного анализа, затраты на передачу данных с Host на Device и обратно не оказывают существенного влияния на общее время, поскольку используется высокочастотная оперативная память и интерфейс PCI Express 3.0 для подключения GPU обеспечивает высокую скорость передачи (до 8 Гбит/с). Однако за счет применения невыгружаемой памяти время передачи данных сокращается примерно в два раза. А при асинхронной передаче удается добиться того, что передача данных осуществляется параллельно с вычислительными операциями на CPU, причем до их завершения. Это позволило, в конечном счете, полностью исключить временные затраты на передачу данных. В целом, передача данных Host — Device занимает больше времени, чем Device — Host, поскольку в первом случае передается матрица, а во втором — вектор.
Проведен анализ времени вычислений при решении прямых задач БКЗ по профилю вдоль скважины длиной 100 м при дискретизации с шагом 0.1м для шести зондов. При этом выполнялась одна факторизация с последующим нахождением 606 решений
Получены оценки времени вычислений при факторизации, нахождении решений СЛАУ и определено общее время на CPU и GPU в зависимости от размера матрицы и количества ее ненулевых элементов. Размер матрицы коэффициентов СЛАУ обусловлен плотностью расчетной сетки, выбор которой зависит от размеров областей неоднородности среды и их относительного контраста УЭС.
Из рис. 4 видно, что при факторизации с увеличением размера матрицы наблюдается ускорение вычислений в случае использования GPU. Так, для матрицы, имеющей порядка 104 элементов, ускорение не наблюдается, а если в матрице порядка 106 элементов, то скорость вычислений увеличивается в три раза. Это обусловлено тем, что
СЛАУ.
104
• GPU
гчзтт
Общее время Решение СЛАУ Факторизация
10-2
104
105
Размер матрицы, эл.
106
Рис. 4. Время вычислений на CPU и GPU при факторизации, нахождении решений СЛАУ и общее время в зависимости от размера матрицы
рассматриваемые матрицы сильно разрежены. Высокое быстродействие на GPU достигается при использовании более плотных матриц, получаемых при численном решении других задач каротажа.
Существенным ограничением при вычислениях на GPU является размер памяти графического устройства. Для хранения матрицы применяется широко используемый формат CSR. В этом формате общий размер памяти для матрицы определяется размером трех векторов. Два из них содержат целочисленные элементы, размер которых задается количеством строк матрицы и количеством ее ненулевых элементов. Третий вектор содержит коэффициенты СЛАУ с одинарной или двойной точностью.
Кроме хранения матрицы необходима память для выполнения факторизации. Так, размер памяти для факторизации матрицы, имеющей порядка 104 элементов, превосходит размер памяти, необходимой для ее хранения, примерно в 20 раз, а для матрицы, имеющей порядка 106 элементов, — в три раза. При этом для вычислений с двойной точностью на GPU в первом случае требуется около 50 МБ, а во втором — 3.5 ГБ. Поэтому при реализации вычислений на GPU необходимо предусматривать достаточный объем памяти для хранения матрицы в формате CSR, а также для ее факторизации.
В отличие от факторизации, нахождение решений СЛАУ на GPU выполняется в несколько раз медленнее, чем на CPU (табл. 1). Так, для матрицы с количеством элементов порядка 104 вычисления на CPU осуществляются в 26 раз быстрее, чем на GPU, а для матрицы с количеством элементов порядка 106 — в четыре раза. Поскольку время передачи данных очень мало, большие временные затраты на поиск решения СЛАУ на GPU обусловлено сильно разреженной матрицей и особенностями вычислений в этом случае. К последним относятся задержки, связанные с обращением к глобальной/общей памяти и синхронизацией потоков. В целом видно, что выигрыш от использования CPU уменьшается с увеличением размера матрицы, но слишком медленно. Исходя из установленной закономерности, вычисления на GPU будут быстрее только для матриц размером более 107 -108 элементов. Все это указывает на то, что нахождение решения СЛАУ эффективнее распараллеливать на CPU.
Нужно отметить, что при решении одной прямой задачи БКЗ (расчеты в заданной точке по профилю скважины для одного зонда), включающей факторизацию и только одно нахождение решения СЛАУ, быстродействие при использовании GPU будет в несколько раз выше. Расчетами установлено, что при нахождении нескольких (от 2 до 13) решений СЛАУ использование GPU будет оправданно. Так, это соответствует решению двух прямых задач БКЗ в двух точках по профилю скважины для шести зондов, что не представляет интереса для практического применения. Но при большом
Таблица 1. Время вычислений на CPU и GPU и общее время в зависимости от размера матрицы и количества ее ненулевых элементов, с
Размер Количество Фактори- Решение Общее
матрицы, ненулевых зация СЛАУ время
тыс. элементов элементов, тыс. элементов CPU GPU CPU GPU CPU GPU
13 119 0.14 0.07 0.67 18.51 1.06 19.03
49 475 0.61 0.32 6.94 69.57 7.47 71.08
192 1897 3.75 1.56 43.54 285.15 50.53 290.70
763 7583 34.64 11.47 339.82 1402.59 386.78 1430.39
Таблица 2. Время вычислений на CPU и GPU при последовательных и параллельных реализациях, а также с использованием гетерогенных вычислений CPU-GPU при решении прямой двумерной задачи БКЗ
Вычислительная операция CPU CPU* GPU CPU-GPU
Факторизация 34.64 8.77 11.47 4.33
Решение СЛАУ 339.82 237.01 1402.59 158.00
количестве решений прямых задач БКЗ (расчеты по профилю вдоль скважины длиной десятки — сотни метров для всех зондов) временные затраты становятся существенными, а скорость вычислений на CPU значительно выше, чем на GPU.
Выполненный анализ показал, что эффективными вычисления будут при факторизации на GPU и нахождении решений СЛАУ на CPU, т. е. при реализации гетерогенных вычислений CPU-GPU. Для разработки параллельного алгоритма применяется библиотека CHOLMOD [29], являющаяся частью вычислительного пакета библиотек SuiteSparse [30] семейства решателей для разреженных матриц, использующих графические процессоры от NVIDIA. Выбор библиотеки CHOLMOD обусловлен тем, что она позволяет выполнять факторизацию и находить решение СЛАУ раздельно, при этом факторизация осуществляется с использованием гетерогенных вычислений CPU-GPU. Поскольку вычисления на CPU и GPU выполняются одновременно, CPU не переходит в длительное ожидание завершения выполнения факторизации на GPU, как это происходит при использовании библиотеки cuSOLVER. Это также позволяет хранить не все данные на GPU, а только необходимые в текущий момент, что обусловливает возможность выполнения вычислений с использованием матриц еще больших размеров, чем рассмотренные в работе.
Данные оценки быстродействия вычислений с использованием разработанных алгоритмов по сравнению с исходным последовательным и параллельным алгоритмами на CPU приведены в табл. 2. Для реализации параллельных вычислений на CPU и GPU применялась библиотека cuSOLVER, а для гетерогенных вычислений CPU-GPU — CHOLMOD. Оценки получены при моделировании данных БКЗ в реалистичной геоэлектрической модели нефтяного пласта сложного геологического строения в случае решения задачи с наибольшим размером матрицы СЛАУ (763 тыс. элементов, 7583 тыс. ненулевых элементов).
Как видно из табл. 2, факторизация с использованием гетерогенных вычислений CPU-GPU выполняется в восемь и два раза быстрее, чем при последовательных и параллельных вычислениях на CPU соответственно, и в 2.6 раза быстрее, чем на GPU. При этом времени на нахождение всех решений СЛАУ с использованием гетерогенных вычислений CPU-GPU затрачивается в 2.2 и 1.5 раза меньше по сравнению с исходным последовательным и параллельным (в табл. 2 отмечено звездочкой) алгоритмами на CPU соответственно и в 8.9 раза меньше по сравнению с распараллеленным на GPU. Экспериментами установлено, что общее время вычислений на CPU-GPU в 1.5-2.3 раза меньше, чем на CPU, и в 8.7 раза меньше, чем на GPU.
Таким образом, полученное быстродействие при гетерогенных вычислениях CPU-GPU демонстрирует высокую эффективность разработанного алгоритма и его программных реализаций, что указывает на возможность их применения для оперативного численного моделирования.
3. Моделирование данных БКЗ на основе гетерогенных вычислений CPU — GPU
С использованием разработанного вычислительного алгоритма и реализованной программы на основе гетерогенных вычислений CPU-GPU проведено моделирование данных БКЗ в геоэлектрической модели сложнопостроенного геологического объекта, содержащего нефтяные пласты. Рассматривается реалистичная геоэлектрическая модель нефтегазового пласта, характерного для Широтного Приобья. Как отмечалось выше, при бурении проницаемых пластов происходит фильтрация в них бурового раствора, что приводит к оттеснению пластового флюида (водонефтяной смеси) и образованию прискважинной зоны. Изменение УЭС в зоне проникновения связано с перераспределением минерального состава пластовой воды и фильтрата бурового раствора. За ней располагается неизмененная часть пласта (рис. 5, а).
a
Рис. 5. Реалистичная геоэлектрическая модель нефтяного пласта сложного геологического строения (а) и синтетические диаграммы кажущегося УЭС зондов БКЗ (б), полученные с использованием гетерогенных вычислений CPU-GPU
Модель среды состоит из двух маломощных (2.0-2.5 м) нефтеводосодержащих пластов и тонкого глинистого слабопроницаемого прослоя (0.5 м) между ними. Верхний нефтенасыщенный пласт имеет значения УЭС в зоне проникновения 20 Ом • м, в неизмененной части пласта — 15 Ом • м, нижний водонасыщенный пласт — 30 и 5 Ом • м соответственно. Толщина прискважинных цилиндрических зон в этих пластах составляет 0.4 и 0.6 м. Пласты отделены друг от друга проводящим глинистым прослоем с УЭС, равным 3 Ом • м.
Анализ диаграмм кажущегося УЭС зондов БКЗ (A0.4M0.1N, A1.0M0.1N, A2.0M0.5N, A4.0M0.5N, A8.0M1.0N), приведенных на рис. 5, б, позволил установить, что они значительно дифференцированы и отражают вертикальную неоднородность геологического разреза, что дает возможность корректно провести его расчленение на пласты. В интервалах проницаемых пластов, характеризующихся проникновением фильтрата бурового раствора и неоднородным строением прискважинной области, наблюдается расхождение диаграмм разноглубинных зондов. Представленные диаграммы соответствуют типичной геоэлектрической ситуации пространственно неоднородных флюидонасы-щенных горных пород. Необходимо отметить, что приведенные результаты полностью идентичны данным, полученным при вычислениях с использованием последовательного алгоритма.
Таким образом, разработанный алгоритм позволит в дальнейшем проводить быстрое численное моделирование сигналов зондов БКЗ и изучать особенности их поведения в сложных геологических условиях, что чрезвычайно важно как при создании новых каротажных технологий и аппаратурных комплексов, так и при оперативной интерпретации практических каротажных данных в масштабе реального времени.
Заключение
Разработан и программно реализован алгоритм моделирования данных электрического каротажа нефтегазовых скважин с использованием высокопроизводительных гетерогенных вычислений CPU-GPU. Алгоритм основан на решении прямой задачи бокового каротажного зондирования в двумерных геоэлектрических моделях среды сеточным методом конечных элементов. Решение системы линейных алгебраических уравнений с разреженной матрицей большой размерности находится с использованием разложения Холецкого.
Программные реализации алгоритма с использованием технологии NVIDIA CUDA и вычислительных библиотек cuSOLVER и CHOLMOD позволяют осуществлять факторизацию и нахождение решения СЛАУ на CPU и GPU как раздельно, так и совместно. Выполнен анализ времени вычислений на CPU и GPU при факторизации, нахождении решений СЛАУ и общего времени вычислений в зависимости от размера матрицы и количества ее ненулевых элементов. Установлено, что применительно к решаемой задаче параллельные вычисления наиболее эффективны при факторизации на GPU и нахождении решений СЛАУ на CPU.
Получены оценки быстродействия вычислений на CPU и GPU, включая высокопроизводительные гетерогенные на CPU-GPU. Установлено, что с использованием гетерогенных вычислений на CPU - GPU удается повысить быстродействие по сравнению с аналогичными расчетами на CPU или GPU. Так, факторизация на CPU-GPU выполняется в 2-8 раз быстрее, чем на CPU, и в 2.6 раза быстрее, чем на GPU. При этом нахождение решений СЛАУ на CPU-GPU в 1.5-2.2раза быстрее, чем на CPU, и в 8.9раза
быстрее, чем на GPU. Экспериментами установлено, что общее время вычислений на CPU-GPU в 1.5-2.3 раза меньше, чем на CPU, и в 8.7 раза меньше, чем на GPU.
Проведено численное моделирование данных БКЗ в сложной геологической среде, содержащей нефтяные пласты, характерные для Западной Сибири. Гетерогенные вычисления, реализованные на CPU -GPU, полученные оценки быстродействия и выполненное моделирование указывают на высокую эффективность разработанного алгоритма для решения широкого круга практических задач.
Список литературы / References
[1] Эпов М.И., Глинских В.Н., Никитенко М.Н., Сухорукова К.В. Комплекс СКЛ для каротажа за один спускоподъем: методики и интерпретация на практических примерах // Пути реализации нефтегазового и рудного потенциала Ханты-Мансийского автономного округа — Югры: Матер. Пятнадцатой науч.-практ. конф. Ханты-Мансийск: ИздатНаукаСервис. 2012. Т. 2. С. 27-36.
Epov, M.I., Glinskikh, V.N., Nikitenko, M.N., Sukhorukova, C.V. SKL complex for logging during one round-trip operation: techniques and practical data interpretation // Ways of Realizing the Petroleum and Ore Potential of the Khanty-Mansi Autonomous Okrug — Yugra: Proc. of the Fifteenth Research-to-practice Conf. Khanty-Mansiysk: IzdatNaukaServis, 2012. Vol. 2. P. 27-36.
[2] Эпов М.И., Глинских В.Н. Быстрое двумерное моделирование высокочастотного электромагнитного поля для задач каротажа // Геол. и геофизика. 2003. Т. 44, № 9. С. 942-952. Epov, M.I., Glinskikh, V.N. Fast two-dimensional simulation of high-frequency electromagnetic field in induction logging // Russian Geology and Geophysics. 2003. Vol. 44, No. 9. P. 904-915.
[3] Эпов М.И., Глинских В.Н. Линеаризация относительных характеристик высокочастотного магнитного поля в двумерных проводящих средах // Геология и геофизика. 2004. Т. 45, № 2. С. 266-274.
Epov, M.I., Glinskikh, V.N. Linearization of relative parameters of a high-frequency magnetic field in two-dimensional conducting media // Russian Geology and Geophysics. 2004. Vol. 45, No. 2. P. 247-257.
[4] Глинских В.Н., Эпов М.И. Локально-нелинейные приближения высокочастотного электромагнитного поля для задач каротажа // Геология и геофизика. 2006. Т. 47, № 8. С. 938-944.
Glinskikh, V.N., Epov, M.I. Locally nonlinear approximation of high-frequency electromagnetic field for logging applications // Russian Geology and Geophysics. 2006. Vol. 47, No. 8. P. 930-936.
[5] Еремин В.Н., Нечаев О.В., Хаберхауэр Ш., Шокина Н.Ю., Шурина Э.П. Параллельная реализация математического моделирования процессов электромагнитного каротажа зондовым комплексом ВИКИЗ // Вычисл. технологии. 2007. Т. 12, № 6. С. 18-33. Eremin, V.N., Nechaev, O.V., Haberkhauer, S., Shokina, N.Yu., Shurina, E.P.
Parallel realization of mathematical modelling of electromagnetic logging processes using VIKIZ probe complex // Comput. Technologies. 2007. Vol. 12, No. 6. P. 18-33. (In Russ.)
[6] Эпов М.И., Шурина Э.П., Архипов Д.А. Параллельные конечно-элементные вычислительные схемы в задачах геоэлектрики // Вычисл. технологии. 2013. Т. 18, № 2. С. 95-112.
Epov, M.I., Shurina, E.P., Arkhipov, D.A. Parallel numerical finite elements schemes for solving of geo-electric problems // Comput. Technologies. 2013. Vol. 18, No. 2. P. 95-112. (In Russ.)
[7] Глинских В.Н., Нестерова Г.В., Эпов М.И. Моделирование и инверсия данных электромагнитного каротажа с использованием петрофизических моделей электропроводности // Геология и геофизика. 2014. Т. 55, № 5-6. С. 1001-1010.
Glinskikh, V.N., Nesterova, G.V. Epov, M.I. Forward modeling and inversion of induction logs from shaly sand reservoirs using petrophysical conductivity models // Russian Geology and Geophysics. 2014. Vol. 55, No. 5-6. P. 793-799.
[8] Эпов М.И., Глинских В.Н., Никитенко М.Н., Сухорукова К.В. Численное моделирование и анализ сигналов электромагнитного каротажа в процессе бурения // Каро-тажник. 2014. Вып. 11(245). С. 29-41.
Epov, M.I., Glinskikh, V.N., Nikitenko, M.N., Sukhorukova, C.V. Numerical simulation and analysis of electromagnetic logging-while-drilling data // Karotazhnik. 2014. Vol. 11(245). P. 29-41. (In Russ.)
[9] Глинский Б.М., Костин В.И., Кучин Н.В., Соловьев C.A., Чеверда В.А. Организация параллельных вычислений для решения уравнения Гельмгольца прямым методом с использованием малоранговой аппроксимации и HSS формата // Вычисл. методы и программирование. 2015. Т. 16. С. 607-616.
Glinskiy, B.M., Kostin, V.I., Kuchin, N.V., Soloyev, S.A., Cheverda, V.A. Aspects of parallel computing to solve Helmholtz equation by a direct solver with low-rank approximation and the HSS format of data storage // Numerical Methods and Programming. 2015. Vol. 16. P. 607-616. (In Russ.)
[10] Костин В.И., Решетова Г.В., Чеверда В.А. Численное моделирование трехмерного акустического каротажа с использованием многопроцессорных вычислительных систем // Матем. моделирование. 2008. Т. 20, № 9. С. 51-66.
Kostin, V.I., Reshetova, G.V., Cheverda, V.A. Parallel numerical simulation of 3D acoustic logging // Matem. Modelirovanie. 2008. Vol. 20, No. 9. P. 51-66. (In Russ.)
[11] Belonosov, M.A., Kostov, C., Reshetova, G.V., Solovev, S.A., Cheverda, V.A.
Parallel numerical simulation of seismic waves propagation with Intel Math Kernel Library // Lecture Notes in Computer Science. Applied Parallel and Scientific Computing. 2013. Vol. 7782. P. 153-167.
[12] Puzyrev, V., Koric, S., Wilkin, St. Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems // Computers and Geosciences. 2016. Vol. 89. P. 79-87.
[13] Глинских В.Н., Эпов М.И., Лабутин И.Б. Моделирование диаграмм электромагнитного каротажа на графических процессорах // Вычисл. технологии. 2008. Т. 13, № 6. С. 50-60.
Glinskikh, V.N., Epov, M.I., Labutin, I.B. Simulation of electromagnetic logs using graphics processing unit // Comput. Technologies. 2008. Vol. 13, No. 6. P. 50-60. (In Russ.)
[14] Глинских В.Н., Буланцева Ю.О. Математическое моделирование данных электромагнитного каротажа с использованием высокопроизводительного сопроцессора Intel Xeon Phi // Вестн. НГУ. Математика, механика, информатика. 2014. Т. 14, вып. 4. С. 11-22. Glinskikh, V.N., Bulantseva, Yu.O. Mathematical simulation of electromagnetic logs using high-performance co-processor Intel Xeon Phi // Vestnik of NGU. Mathematics, Mechanics, Informatics. 2014. Vol. 14, iss. 4. P. 11-22. (In Russ.)
[15] Глинских В.Н., Горбатенко В.А. Инверсия данных электромагнитного каротажа на графических ускорителях // Вычисл. технологии. 2015. Т. 20, № 1. С. 25-37. Glinskikh, V.N., Gorbatenko, V.A. Electromagnetic logging data inversion on GPU // Comput. Technologies. 2015. Vol. 20, No. 1. P. 25-37. (In Russ.)
[16] Labutin, I.B., Surodina, I.V. Algorithm for sparse approximate inverse preconditioners in the conjugate gradient method // Reliable Computing. 2013. Vol. 19. P. 120-126.
[17] Dudaev, A.R., Nechaev, O.V., Glinskikh, V.N. High-performance computing on GPU for 2D electric logging problem based on the finite element methods // Proc. of the 8th Intern. Siberian Early Career GeoSci. Conf. Novosibirsk: NSU, 2016. P. 390-391.
[18] Суродина И.В., Нестерова Г.В. Моделирование показаний зондов ВИКИЗ и БКЗ на графических процессорах. Петрофизика сложных коллекторов: проблемы и перспективы 2015: Сб. статей. М.: ЕАГЕ Геомодель, 2015. С. 85-94.
Surodina, I.V., Nesterova, G.V. Simulation of VIKIZ and BKZ logs on GPUs. Petrophysics of complex reservoirs: challenges and opportunities 2015: Collection of Articles. Мoscow: EAGE Geomodel, 2015. P. 85-94. (In Russ.)
[19] Epov, M.I., Sukhorukova, C.V., Glinskikh, V.N., Nikitenko, M.N., Nechaev, O.V., Surodina, I.V. Effective electromagnetic log data interpretation in realistic reservoir models // Open J. of Geology. 2013. Vol. 3, iss. 2B. P. 81-86.
[20] Davis, T.A., Rajamanickam, S., Sid-Lakhdar, W.M. A survey of direct methods for sparse linear systems // Acta Numerica. 2016. Vol. 25. P. 383-566.
[21] Pierce, D.J., Hung, Y., Liu, C.-C., Tsai, Y.-H., Wang, W., Yu, D. Sparse multifrontal performance gains via NVIDIA GPU // Work. on GPU Supercomputing. Taipei: National Taiwan University, 2009. Available at:
http://cqse.ntu.edu.tw/cqse/download_file/DPierce_20090116.pdf (accessed 30.09.16).
[22] Lacoste, X., Ramet, P., Faverge, M., Ichitaro, Y., Dongarra, J. Sparse direct solvers with accelerators over DAG runtimes: Research Report RR-7972, INRIA. Bordeaux, France, 2012. P. 11. Available at: https://hal.inria.fr/hal-00700066v2/document (accessed 30.09.16).
[23] Mittal, S., Vetter, J.S. A survey of CPU -GPU heterogeneous computing techniques // ACM Computing Surveys. 2015. Vol. 47, iss. 4. P. 69:1-69:35.
[24] Sao, P., Vuduc, R., Li, X.S. A distributed CPU-GPU sparse direct solver // Proc. Euro-Par 2014 Parallel Processing / F. Silva et al. (Eds). Lecture Notes in Computer Science. 2014. Vol. 8632. P. 487-498.
[25] Дахнов В.Н. Электрические и магнитные методы исследования скважин. М.: Недра, 1981. 344 с.
Dakhnov, V.N. Electrical and magnetic borehole survey techniques. Moscow: Nedra, 1981. 344 p. (In Russ.)
[26] Шайдуров В.В. Многосеточные методы конечных элементов. М.: Наука, 1989. 288 с. Shaidurov, V.V. Multigrid finite element methods. Moscow: Nauka, 1989. 288 p. (In Russ.)
[27] CUDA C Programming Guide. Design Guide // NVIDIA CUDA. 2016. Available at: http://docs.nvidia.com/cuda/pdf/CUDA_C_Programming_Guide.pdf (accessed 30.09.16).
[28] cuSOLVER // NVIDIA CUDA Toolkit Documentation. 2016. Available at: http://docs.nvidia.com/cuda/cusolver (accessed 30.09.16).
[29] Chen, Y., Davis, T.A., Hager, W.W., Rajamanickam, S. Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate // ACM Trans. Math. Softw. 2008. Vol. 35(3). P. 22:1-22:14.
[30] SuiteSparse: A suite of sparse matrix software v.4.5.3 / Davis, T.A. SuiteSparse. 2016. Available at: http://www.suitesparse.com (accessed 30.09.16).
Поступила в редакцию 20 октября 2016 г., с доработки — 6 марта 2017 г.
High-performance CPU —GPU heterogeneous computing in resistivity logging of oil and gas wells
Glinskikh, Viacheslav N.1'2'*, Dudaev, Alexander R.1,2, Nechaev, Oleg V.1
1Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia 2 Novosibirsk State University, Novosibirsk, 630090, Russia
* Corresponding author: Glinskikh, Viacheslav N., e-mail: GlinskikhVN@ipgg.sbras.ru
The work is concerned with the development of numerical algorithms for solving forward problems of borehole geoelectrics by applying high performance GPU computing.
Studying new types of complex hydrocarbon reservoirs demands the improvement of logging tools and software for data processing. The improvement of processing accuracy is feasible with the use of numerical solutions in full mathematical formulations. Computational problems in resistivity logging are known to be resource-intensive and not very effective when used in practice. This research considers one of the ways to speed up calculations, namely through the application of NVIDIA graphics processors.
We have developed and implemented into software an algorithm for the simulation of resistivity logging data from oil and gas wells, by making use of high-performance CPU -GPU heterogeneous computations. The numerical solution of the 2D forward problem is based on the finite-element method and the Cholesky decomposition for solving a system of linear equations. The software implementations of the algorithm are made by means of NVIDIA CUDA technology and computing libraries making it possible to decompose the equation system and find its solution on CPU and GPU. The analysis of computing time as a function of the order of matrix and number of non-zero elements has shown that in the case at hand the computations are the most effective when decomposing on GPU and finding a solution on CPU. We have estimated the operational speed of CPU and GPU computations, including high-performance heterogeneous CPU-GPU ones. It is found that heterogeneous CPU-GPU computations enable speeding up in comparison with similar CPU or GPU calculations. Using the developed algorithm, we have simulated resistivity data in realistic models. The results of our investigation point to a high efficiency of the algorithm in respect to dealing with a wide variety of practical problems.
Keywords: graphics processing units, CPU-GPU heterogeneous computing, parallel algorithm, finite-element method, 2D forward problem, resistivity logging.
Received 20 October 2016 Received in revised form 6 March 2017
© ICT SB RAS. 2017
