Для задачи гидродинамического моделирования нефтегазовых месторождений характерна высокая вычислительная ресурсоемкость. Настоящая работа нацелена на ускорение соответствующих расчетов посредством использования гибридных вычислительных систем с графическими процессорами, главным образом, на ускорение решения линейных систем, возникающих при численном моделировании многофазной фильтрации. Исследовано влияние формата хранения матриц на время выполнения на графическом процессоре базовых операций предобусловленного метода бисопряжённых градиентов со стабилизацией. Продемонстрирована производительность решения тестовых разреженных систем при использовании различных предобуславливателей и подходов к их распараллеливанию.
2. Алгоритмы решения СЛАУ и их программная реализация Матрицы СЛАУ, возникающие в ходе численного решения уравнений многофазной фильтрации потоков углеводородов в пористой среде, являются сильно разрежёнными и плохо обусловленными, поэтому для их решения принято использовать итерационные методы подпространства Крылова, например, метод бисопряжённых градиентов со стабилизацией (BiCGStab) с различными предобуславливателями. Из-за высокой жесткости системы найти решение без предобуславливателя или с простейшими вариантами (например, Якоби) не удается, поскольку при большом числе итераций накопление вычислительной ошибки приводит к значительному искажению приближенного решения [2].
Одним из наиболее часто используемых предобуславливателей, применяемых при
численном решении уравнений многофазной фильтрации углеводородов, является неполное
LUразложение без заполнения (с нулевым заполнением) –– ILU(0) [
Алгоритмы построения предобуславливателя ILU(0) и решения сопутствующих систем с
нижне- и верхнетреугольными матрицами в классической постановке обладают
незначительным ресурсом параллелизма для СЛАУ с сильно разреженными матрицами. Существует
несколько различных подходов к распараллеливанию предобуславливателей класса ILU. Одним
из таких подходов является разделение матрицы на уровни, содержащие строки матрицы,
которые могут обрабатываться параллельно (level scheduling) [
В нашем решателе заложены различные подходы к распараллеливанию ILU(0).
1) На CPU используется блочный метод Якоби, обеспечивающий возможность
параллельного исполнения операций построения и применения неполного
LUразложения на множестве процессорных ядер.
2) На одном GPU средствами библиотеки cuSPARSE реализуется разделение строк
матрицы на уровни и их дальнейшая параллельная поуровневая обработка.
3) На одном или нескольких GPU (в рамках CPR) предлагается использовать
комбинированный подход: процедура разделения на уровни применяется к
блочнодиагональной части исходной матрицы с произвольно заданным количеством
блоков (в нашей работе от 64 до 256) [
CPU-версия решателя является многопоточной (OpenMP) и базируется на функциях математической библиотеки Intel Math Kernel Library (MKL). Для хранения матриц используется формат CSR, т.к. другие форматы в функциях MKL, связанных с ILU(0), не поддерживаются. Реализован итерационный метод BiCGStab с предобуславливателем ILU(0), распараллеливание которого проводится как отмечено выше в пункте 1. GPU-версия решателя разработана средствами CUDA, OpenMP, а также базируется на функциях математических библиотек NVIDIA cuBLAS и cuSPARSE из состава CUDA Toolkit и свободной для некоммерческого использования библиотеки AmgX. Реализован итерационный метод BiCGStab с предобуславливателями ILU(0) и CPR. Причем распараллеливание ILU(0) осуществляется в соответствии с приведенными выше пунктами 2-4 в зависимости от конфигурации запуска. Операции построения и применения AMG на данный момент выполняются только на одном GPU, т.к. параллельная версия классического алгоритма алгебраического многосеточного метода для систем с несколькими GPU, реализованная в библиотеке AmgX, может быть задействована только из MPI-программ. Для хранения разреженных матриц в памяти поддерживается как универсальный формат CSR (Compressed Sparse Row), так и формат BSR (Block compressed Sparse Row), предназначенный для хранения разреженных матриц с блочной структурой. В целях обеспечения максимальной производительности GPU-версии линейного решателя далее в разделе 3.1 исследуется влияние формата хранения матриц на время выполнения на GPU базовых вычислительных операций. 3. Экспериментальная часть
В таблице 1 рассмотрены характеристики тестовых разреженных матриц СЛАУ, полученных при гидродинамическом моделировании реальных нефтегазовых месторождений. Отметим, что тестовые матрицы получены при решении уравнений трехфазной фильтрации, дискретизированных по времени с использованием полностью неявной разностной схемы, и для них характерна группировка ненулевых элементов в блоки размером 3x3, что учитывается при использовании формата хранения BSR. Особенности структуры подобных матриц более подробно рассмотрены в [14]. Данные матрицы будут использованы далее для оценки различных аспектов эффективности решения СЛАУ на CPU и GPU.
Таблица 1. Характеристики тестовых матриц 3.1 Исследование влияния формата хранения матриц на время выполнения на GPU базовых операций предобусловленного метода BiCGStab
К наиболее трудоемким (базовым) операциям метода BiCGStab с предобуславливателем ILU(0) относятся следующие: • построение предобуславливателя ILU(0); • решение треугольных систем; • умножение матрицы на вектор.
Ранее в нашем решателе как в CPU-, так и в GPU-версиях для хранения разреженных матриц использовался популярный формат CSR, который поддерживается библиотеками MKL и cuSPARSE во всех указанных выше операциях. С выходом CUDA Toolkit версии 6.0 в библиотеке cuSPARSE появилась поддержка формата BSR в функциях, связанных с использованием неполного LU-разложения, причем для операции умножения разреженной матрицы на вектор соответствующая поддержка имелась и ранее. В связи с этим в GPU-версии решателя в качестве одного из основных форматов хранения матриц был поддержан формат BSR. В данном разделе приводятся результаты исследования влияния используемого формата хранения матриц на производительность выполнения на GPU упомянутых выше базовых операций.
Необходимо отметить, что на данный момент библиотека cuSPARSE содержит две версии функций, реализующих построение ILU(0) и решение треугольных систем, для разреженных матриц в формате CSR. Например, *csrilu0 соответствует первой версии функции построения ILU(0), а *csrilu02 – второй версии. С учетом того, что для формата BSR первые версии функций отсутствуют, далее в сравнении используются только вторые версии функций и для формата CSR.
В таблице 2 приведены времена выполнения базовых операций при работе с различными форматами хранении матриц. Видно, что на операциях, связанных с ILU(0) наибольшая производительность достигается при использовании формата BSR: время построения ILU(0) меньше в 2,3-2,6 раза, а время решения треугольных систем меньше в 1,2-1,5 раза. В то же время операция умножения разреженной матрицы на вектор выполняется быстрее при использовании формата CSR, соответствующее ускорение составляет от 2 до 3 раз.
Таблица 2. Время выполнения на GPU базовых операций метода BiCGStab с
предобуславливателем ILU(0) 3.2 Оценка эффективности решения СЛАУ на CPU и GPU методом BiCGStab с предобуславливателем ILU(0)
В данном разделе приводится сравнительная оценка эффективности решения тестовых СЛАУ на 1-16 ядрах CPU и 1-2 GPU. Отметим, что распараллеливание операций построения и применения ILU(0) для CPU производится с помощью блочного метода Якоби: матрица, на основе которой строится ILU(0), разбивается на диагональные блоки (по числу задействованных процессорных ядер), далее обрабатываемые независимо. Аналогично производится разделение работы между несколькими GPU. Данный подход приводит к утрате части ненулевых элементов в матрице предобуславливателя с увеличением количества блоков, что негативно влияет на сходимость итерационного метода. Это можно заметить в таблице 3: количество итераций при расчете на 16 ядрах CPU больше в 1,04-3,24 раза относительно последовательного расчета. Таблица 3. Количество итераций метода BiCGStab с блочно-диагональным предобуславливателем
ILU(0) при решении СЛАУ на CPU и GPU Время решения СЛАУ, включающее в себя время на построение предобуславливателя и итерации метода BiCGStab, на 1-16 ядрах CPU и 1-2 GPU приведено в таблице 4. Максимальное ускорение параллельных вычислений на CPU составляет от 2,14 до 6,7 раз в зависимости от матрицы СЛАУ, причем минимальные ускорения наблюдаются на матрицах, для которых характерен рост числа итераций с увеличением количества задействованных ядер CPU. Ускорение при расчете на 2 GPU относительно 1 GPU составляет 1,61-1,73 раза.
Таблица 4. Время решения СЛАУ на CPU и GPU методом BiCGStab с блочно-диагональным предобуславливателем ILU(0) Целесообразность использования GPU подтверждается снижением времени решения СЛАУ на 2 GPU в 2,23-4,03 раза относительного минимального времени расчета на 2 CPU. Рис. 1. Ускорение решения СЛАУ на CPU и GPU методом BiCGStab с блочно-диагональным предобуславливателем ILU(0) Ускорение, которое можно получить при параллельном решении различных тестовых СЛАУ методом BiCGStab с блочно-диагональным предобуславливателем ILU(0), относительно времени последовательного расчета на одном ядре CPU продемонстрировано на рисунке 1. 3.3 Оценка эффективности решения СЛАУ на GPU методом BiCGStab с предобуславливателем CPR
Таблица 5 иллюстрирует улучшение сходимости итерационного метода при использовании предобуславливателя CPR: например, при решении СЛАУ с матрицей krrv число итераций может быть снижено до 7 раз. При этом возможно без ухудшения сходимости применение в рамках CPR комбинированного подхода к параллельному построению неполного LU-разложения без заполнения (ILU(0)*), который предполагает построение предобуславливателя в соответствии с алгоритмом разделения на уровни на базе блочно-диагональной части исходной матрицы с заданным количеством блоков. Таблица 5. Количество итераций метода BiCGStab с различными предобуславливателями при решении СЛАУ на GPU ILU(0)
Количество итераций
CPR (AMG + ILU(0)) 1 GPU 2 GPU 1 GPU 2 GPU 1 GPU В таблице 6 приведены времена решения СЛАУ на 1-2 GPU при использовании различных предобуславливателей: ILU(0), CPR и CPR*, в рамках которого реализуется комбинированный подход при работе с ILU(0) для извлечения дополнительного параллелизма. Видно, что снижение количества итераций, отмеченное выше при использовании CPR, преимущественно приводит к снижению времени, затрачиваемого на 1 GPU на решение СЛАУ (построение предобуславливателя и итерационный процесс). Исключением является матрица imsh, для которой использование CPR приводит к замедлению решения СЛАУ, т.к. построение предобуславливателя CPR является достаточно трудоемкой операцией, в то время как для обеспечения хорошей сходимости метода (1,5 итерации) достаточно предобуславливателя ILU(0).
Таблица 6. Время решения СЛАУ на GPU методом BiCGStab с различными предобуславливателями ILU(0)
Время решения СЛАУ, с
CPR (AMG + ILU(0)) 1 GPU 2 GPU 1 GPU
2 GPU Ускорение решения СЛАУ с предобуславливателями CPR и CPR* на 2 GPU относительно 1 GPU составляет от 1,2 до 1,34 раза. Относительно низкая масштабируемость связана с тем, что значительную часть от общего времени решения СЛАУ (не менее 24%) при расчете на 1 GPU занимает классический алгебраический многосеточный метод, реализация которого в библиотеке AmgX не позволяет загрузить несколько GPU из многопоточной программы. Применение комбинированного подхода к распараллеливанию ILU(0) в рамках CPR позволяет снизить время решения СЛАУ на 3-16% для одного GPU и 3-12% для двух GPU. 4. Заключение
Продемонстрированная в работе производительность решения СЛАУ, возникающих в ходе численного решения уравнений многофазной фильтрации потоков углеводородов в пористой среде, позволяет говорить о возможности применения GPU NVIDIA в рамках решения задачи гидродинамического моделирования нефтегазовых месторождений. В целях повышения масштабируемости решателя на системах с несколькими GPU в перспективе планируется разработка MPI-версии, которая обеспечит возможность параллельной работы алгебраического многосеточного метода на нескольких GPU. Литература Development of parallel linear solver for reservoir simulation on hybrid computing systems with GPUs Arthur Yuldashev, Ratmir Gubaidullin and Nikita Repin Keywords: graphics processors, sparse linear systems, multi-and many-core systems, parallel computing, reservoir simulation Problem of hydrodynamic modeling of oil and gas fields is characterized by high consumption of computing resources. The present work aims to accelerate reservoir simulations through the use of hybrid computing systems with graphics processors, mainly to accelerate linear solver for numerical simulations of multiphase filtration.