При распараллеливании решения трехмерных краевых задач на квазиструктурированных сетках методом декомпозиции расчетной области на подобласти, сопрягаемые без наложения, наиболее трудоемкой вычислительной процедурой является решение краевых подзадач в подобластях. Использование параллелепипедальных квазиструктурированных сеток дает возможность применить для этих целей быстросходящиеся методы переменных направлений. Распараллеливание итерационного процесса по подобластям проводится на CPU в системе MPI, а для решения подзадач в настоящей работе предлагается использовать графические ускорители GPU. Проводятся экспериментальные исследования применения графических ускорителей при решении подзадач методом Писмана - Рэчфорда. Даются экспериментальные оценки ускорения распараллеливания в гибридной вычислительной среде CPU + GPU по сравнению с расчетами только на CPU.
Исследование эффективности применения графических ускорителей при распараллеливании решения трехмерных краевых задач на квазиструктурированных сетках* нологии его реализации в системе CUDA [8] и даются результаты численных экспериментов, показывающие значительное (более 60 раз) ускорение вычислений по сравнению с расчетами только на CPU. 2. Постановка задачи и основы алгоритма ее решения
Пусть в замкнутой трехмерной области G G с границей требуется решить краевую задачу
u g1 , l u g2 . (1) Здесь u u(T ) – искомая функция, g1 g1(T ) , g2 g2 (T ) – заданные функции ( T (x, y, z) – текущая точка, где x, y, z – декартовы координаты), – оператор Лапласа, l – оператор граничных условий. Рассматриваются граничные условия Дирихле, Неймана, а также смешанные краевые условия. Предполагается, что граница и функции g1 , g2 таковы, что существует единственное решение задачи (1), обеспечивающее гладкость, достаточную для проведения дальнейших рассуждений.
Построим в расчетной области G структурированную равномерную макросетку H с шагами, намного превышающими максимальный шаг результирующей сетки, на которой ищется решение исходной задачи. Тем самым мы проведем декомпозицию G на непересекающиеся подобласти Gm , m 1, M ,где M – известное целое число. Граница сопряжения подобластей (интерфейс) , в свою очередь, разбивается на грани f , ребра e и макроузлы m , являющиеся узлами макросетки H так, что f e m . Для дальнейшего удобно ввести объединение e,m e m и область G0 G \ .
Построим в подобластях Gm структурированные равномерные подсетки h,m . Объединение этих подсеток составляет квазиструктурированную сетку
M h h,m .
m1 Исходную краевую задачу (1) переформулируем следующим образом: в замкнутой области G требуется найти решение уравнения Пуанкаре – Стеклова
Fv(T ) 0, T f , (2) совместно с решением краевых задач
u(T ) g1(T ), T e,m (3) u(T ) g1(T ) , l u g2 , u v, T G0 (4) относительно функций u и v . Здесь оператор F определяется как
uv() uv() Fv , (5)
n f n f где v – след функции u на (в том числе на f ).
Решение задачи (2) – (4) будем проводить методом итераций по подобластям, состоящим из следующих этапов.
1. Задается начальное приближение v(f0) на гранях f . 2. Из уравнения (3) находятся значения функции ue(,nm) ve(,nm) ( n 0,1, – номер итерации) на ребрах и в макроузлах e,m . 3. Из решения краевой задачи (4) с граничными условиями Дирихле v v(n) v(fn) ve(,nm) находятся значения искомой функции на n - ом приближении в подобластях. 4. Рассчитываются производные, входящие в выражение (5). 5. Делается очередной (n 1) -ый шаг по решению уравнения Пуанкаре – Стеклова (2) и находятся значения функции v(fn1) на гранях. 6. Если сходимость итерационного процесса достигнута, то находится окончательное решение на ребрах, в макроузлах и в подобластях, если же это не так, то повторяется пункты 2 – 5.
На квазиструктурированной сетке h краевая задача (4) методом конечных разностей, конечных элементов или конечных объемов заменяется приближенной задачей h uh g1 , lh uh g 2 , uh vh , где uh , vh – приближенные значения функций u , v , а h , lh – аппроксимации оператора Лапласа и оператора граничных условий. Ее решение сводится к решению дискретных подзадач на подсетках h,m , что может быть выполнено параллельно. Тем самым мы реализуем на каждой итерации по подобластям этап 3, описанного выше алгоритма.
Для решения подзадач на интерфейсе введем на гранях f сетку f , на ребрах e – сетку e . Для единообразия сетку из макроузлов обозначим как m m . В объединении e,m будем рассматривать сетку e,m e m .
Для нахождения функции v f заменим в (5) производные приближенными разностными соотношениями и потребуем, чтобы разность приближенных производных в узлах сетки f обращалась в нуль, что дает систему линейных алгебраических уравнений (подробнее см. [3])
Av f b 0 , где A – квадратная матрица, а b , v f – векторы.
Элементы матрицы A и вектора b не известны. Известно лишь действие оператора F , h аппроксимирующего оператор F , на какую-либо функцию v~f , заданную в узлах f . Это действие определяется формулой
Fhv~f Av~f b .
Для вычисления компонент вектора b дадим v~f пробное значение v~f(0) 0 , решим соответствующую краевую задачу, а затем, вычислив Fhv~f(0) , получим искомый вектор.
Решение системы линейных алгебраических уравнений (7) будем проводить каким –либо итерационным методом в подпространствах Крылова [7]. Замечательным свойством данных методов, которое мы используем, является то, что они не требуют знания элементов матрицы A , а требуют лишь действия A на некий вектор p , что согласно (8) может быть вычислено как
Ap Fh p b . На каждом шаге крыловского итерационного процесса необходимо решать краевые подзадачи в подобластях, поэтому его называют итерационным процессом по подобластям. Положительным свойством трудоемкой процедуры итераций по подобластям является то, что она допускает параллельную реализацию.
Решение подзадач на ребрах и в макроузлах составляет этап 2, описанного выше алгоритма и проводится на каждой n -ой итерации по подобластям. К моменту проведения данных вычислений значения функции v(fn) на гранях считаются известными. (6) (7) (8) и, соответственно, векторы, входящие в (9), как (9) (10) Здесь блоки имеют следующие размерности: B11 : Ne Ne , B12 : Ne Nm , B21 : Nm Ne , B22 : Nm Nm , где N e – число реберных узлов e , а Nm – число макроузлов m, .
Решение системы (9) будем проводить при помощи итерационного процесса, каждый шаг которого состоит из двух полушагов
B11v~e 12 B12v~m fe B21v~e 12 B22v~m 1 fm , v~ 1 v~e 12 , e где 0,1, – номер итерации.
На первом из них вычисляются значения функции на ребрах при фиксированных значениях в макроузлах, а на втором – корректируются значения в макроузлах при фиксированных значениях на ребрах. На ребрах мы имеем «одномерные» сеточные уравнения, так как известные значения функции v(fn) на гранях, входящие в аппроксимацию уравнения (3), исключаются и переносятся в правую часть fe .
Второй полушаг по сути дела состоит в пересчете значений в отдельных макроузлах, которые также не связаны друг с другом, а связаны с полученными на первом полушаге значениями на ребрах.
Рассмотрим решение подзадач в подобластях. Запишем подзадачу вида (6) в m - ой подобласти в матричном виде
Du g , (11) где D – квадратная матица, u – искомый вектор (индексы h, m мы опускаем), g – известный вектор. Рассмотрим аналог итерационного неявного метода Писмана–Речфорда для трехмерного случая. При решении системы (11) он может быть представлен как
un1/ 2 un1 z Dxyun1/ 2 Dzun1 g, un un1/ 2 z Dxyun1/ 2 Dzun g, (12) где n 1, 2, – номера итераций, z – итерационный числовой параметр, а Dxy и Dz – пятидиагональная и трехдиагональная матрицы такие, что D Dxy Dz (подробнее см. [7]). Отсюда видно, что на полуцелом шаге решаются «двумерные» системы, а на целом – «одномерные». Решение «двумерных» систем также проводится методом Писмана Рэчфорда вида (12), в котором матрица Dxy представляется в виде Dxy Dx Dy , где D , D y – трехдиагональные матx рицы. Решение всех «одномерных» систем, к которым приводит реализация алгоритма вида (12), осуществляется методом прогонки. На сетке e,m с привлечением узлов сетки f аппроксимируем исходное уравнение (3). Получим систему сеточных уравнений с матрицей B относительно функции v~ , определенной на ребрах и в макроузлах. Матрицу B можно представить в блочном виде
Bv~ f
B11 B B21 v~ ~v~e , f ffme .
vm 3. Технологии распараллеливания
При распараллеливании рассматриваемой задачи применяется две парадигмы – MPIраспараллеливание и CUDA-распараллеливание, реализуемые соответственно на вычислительном кластере и вычислительных ускорителях GPU. 3.1 MPI-распараллеливание
MPI-распараллеливанию, в первую очередь, подлежит итерационный процесс по подобластям, который занимает подавляющую часть времени решения всей задачи, и, во вторую очередь, внутренний итерационный процесс (10) поиска решений на ребрах и в макроузлах. Эти итерационные процессы обладают внутренним естественным параллелизмом и не требуют дополнительных вычислительных затрат. Эффективность их распараллеливания целиком и полностью зависит от технологии проведения расчетов, от организации обменов между процессорами вычислительной сети.
Важным звеном в технологической цепи решений по достижению эффективности распараллеливания является отображение сеточных данных на вычислительную сеть. Обычно принятый способ отображения: «одна подобласть – один процессор» в данном случае не эффективен, так как подобласти могут содержать различное число узлов, в которых вычисляются значения искомой функции (в дальнейшем – счетных узлов), что приводит к разбалансировке загрузки процессоров. Поэтому подобласти группируются в объединения, содержащие приблизительно одинаковое число счетных узлов, и устанавливается отображение: «одно объединение – один процессор». Передача информации между подобластями одного объединения обходится без межпроцессорных обменов, а между подобластями различных объединений требуется их проведение. Инициализация и осуществление межпроцессорных обменов являются самыми медленными операциями процесса расчета. В связи с этим в основу технологии MPIраспараллеливания кладутся асинхронные операции, позволяющие хотя бы частично проводить обмены на фоне выполнения необходимых арифметических и логических опреаций. 3.2 CUDA-распараллеливание
Решение подзадач в каждой подобласти осуществляется при помощи трехмерного аналога метода Писмана – Речфорда, который распараллеливается на вычислительных ускорителях GPU.
CUDA – это архитектура параллельных вычислений от NVIDIA [8], позволяющая существенно увеличить вычислительную производительность благодаря использованию графических процессоров GPU (Graphics Processing Unit). GPU состоит из процессорных ядер, которые в терминологии NVIDIA называются Streaming Multiprocessor (SM), каждое из которых выполняет сотни программных нитей. Выполнение нитей на GPU происходит блоками, каждый из которых может состоять от 1 до 1024 нитей. Объединение исполняемых блоков называется гридом (Grid). В CUDA блок можно представить в виде одно-, двух- или трехмерного массива, где индекс массива – это индекс исполняемой нити. Множество блоков в гриде представляется в виде одного двумерного массива.
Программу, реализующую трехмерный аналог метода Писмана – Речфорда, можно условно разделить на три части: 1) Выделение памяти на GPU и копирование данных из памяти CPU. 2) Запуск итерационного процесса на GPU. 3) Копирование результата в память CPU и удаление выделенной памяти на GPU. Копирование и выделение данных на GPU осуществляется средствами CUDA. Данные, являющиеся константами в итерационном процессе, копируются в константную память устройства до проведения итераций. Блок – схема проведения итерационного процесса показана на рис. 1. Рис. 1.Блок схема итерационного процесса Рассматриваемый алгоритм реализован в виде функции языка C. На ее вход подается информация о краевых условиях и правой части исходной краевой задачи, параметры подсетки, на которой решается рассматриваемая в данной подобласти краевая подзадача, счетные параметры:
Работу данной функции можно условно разделить на две части: 1) подготовка и 2) расчет. В первой части вычисляются коэффициенты для метода прогонки; выделяется память на GPU, в которую копируются коэффициенты прогонки, краевые значения и значения правых частей уравнения Пуассона. Во второй части реализован двухуровневый итерационный процесс в виде двухуровневого цикла. Во вложенном цикле на GPU запускаются два расчета одномерных задач вдоль осей OX и OY и расчет невязки двумерной задачи, по окончании которых проверяется условие на выход из цикла. В основном цикле на GPU запускается расчет одномерных задач вдоль оси OZ и расчет невязки трехмерной задачи. 4. Численные эксперименты
Цель проводимых численных экспериментов – исследование эффективности применения графических ускорителей при распараллеливании решения трехмерных краевых задач на квазиструктурированных сетках. Критерием эффективности служила величина отношения времени расчетов с использованием только CPU ко времени расчетов с использованием CPU и GPU. Рассматривалась следующая модельная задача
u 0, u 1 в единичном кубе R . Проводилась декомпозиция расчетной области при помощи равномерной макросетки
1 1 1 H X I IH x ,YJ JH y , Z K KH z , I 0, N x , J 0, N y , K 0, N z , H x , H y , H z , N x N y N z где N x , N y , N z – заданные целые числа.
В замкнутых макроэлементах Rm RI,J,K строились равномерные параллелепипедальные подсетки h,m xim X I imhx,m , y jm YJ jmhy,m , zkm ZK kmhz,m im 0, nx,m , jm 0, ny,m , km 0, nz,m с шагами hx,m
X I 1 X I , nx,m hy,m YJ 1 YJ , hz,m Z K 1 Z K .
ny,m nz,m В качестве оператора h выбирались обычные семиточечные разностные операторы Лапласа. Итерации по подобластям проводились при помощи метода сопряженных градиентов, легко поддающегося распараллеливанию.
Вычисления проводились на кластере НКС-30Т Сибирского суперкомпьютерного центра (ССКЦ СО РАН, г. Новосибирск), который имеет следующие характеристики одного вычислительного узла: • CPU:
2 CPU Xeon X5670 6 х 2.93 GHz • GPU: 3 NVIDIA Tesla M 2090 6Gb на архитектуре Fermi 512 CUDA Cores Проведенные серии экспериментов тематически можно представить в виде трех разделов, приведенных ниже. 4.1 Исследование эффективности в зависимости от числа ядер CPU В таблице 1 приведено отношение времени, затраченного на вычисление решения только на CPU ядрах к времени, затраченному на вычисление решения с использованием GPU. Макросетка имела параметры N x , N y 4, N z 3 , а параметры согласованных подсеток приведены в таблице 1. При запуске расчета на GPU устройстве размер блока брался равным 8x8, что позволило обеспечить эффективное использование параметров GPU.
Из полученных результатов можно сделать следующие выводы: 1) ускорение, полученное с использованием GPU, с точностью до погрешностей измерения не зависит от количества ядер CPU, 2) рост ускорения с ростом числа узлов в сетке связан с возрастающим объемом вычислений, приходящихся на одну нить GPU устройства.
Таблица 1. Ускорение в зависимости от количества ядер CPU Ядра CPU 1 3 6 12 5 4 4 4 20 19 18 19
Сетка 39 37 37 38 33 34 34 36 45 47 47 47 4.2 Исследование эффективности на согласованных и несогласованных квазиструктурированных сетках
В проводимых экспериментах макросетка имела параметры N x , N y 4, N z 3 , а размеры подсеток по одному направлению приведены на рис.2. Несогласованная сетка была такой, что любые две соседние подсетки имели различное, отличающееся в два раза по одному направлению, число узлов. На рис. 2 дано число узлов густых подсеток в несогласованных сетках. Рис. 2. Ускорение, полученное на различных сетках на 3 ядрах CPU Из приведенного рисунка следует, что эффективность применения GPU незначительно отличается для согласованных и несогласованных сеток. Количество блоков, вычисляемое на GPU, зависит от размера сетки. Если это количество не кратно количеству блоков, которое GPU устройство может вычислять одновременно, то это ведет к снижению эффективности вычисления. Данное свойство объясняет нелинейность ускорения, изображенного на рис. 2. 4.3 Исследование эффективности в зависимости от размеров макросетки На рис. 3 показана зависимость полученного ускорения в зависимости от размера сетки в подобласти, при различных размерах макросетки. Различие в ускорениях объясняется тем, что, с ростом количества узлов макросетки, растет время на обмен и обработку данных между подобластями на одном ядре CPU.
Рис. 3. Ускорение при разном количестве узлов макросетки 5. Заключение
В работе исследована эффективность применения графических ускорителей при распараллеливании решения трехмерных краевых задач на квазиструктурированных сетках. Распараллеливание осуществляется методом декомпозиции расчетной области на подобласти, сопрягаемые без наложения, основанном на прямой конечно-разностной аппроксимации уравнения Пуанкаре – Стеклова на интерфейсе. Возникающий при этом итерационный процесс по подобластям распараллеливается на CPU. Для решения подзадач в подобластях, составляющего наибольшую часть времени решения всей задачи, применяется метод Писмана – Рэчфорда, обладающий быстрой сходимостью. Его распараллеливание осуществляется на GPU в системе CUDA. Показано, что применение графических ускорителей значительно (более 60 раз) сокращает время решения задачи по сравнению с расчетами только на CPU. Литература 1. 2. 3.
Василевский Ю.В., Ольшанский М.А. Краткий курс по многосеточным методам и методам декомпозиции области. М.: МГУ. 2007.
Quarteroni A., Valli A. Domain Decomposition Methods for Partial Differential Equations // Oxford: Clarendon Press. 1999. Свешников В.М. Построение прямых и итерационных методов декомпозиции // СибЖИМ. 2009. Т.12, № 3(39). С. 99 – 109. Свешников В.М., Беляев Д.О. Построение квазиструктурированных локальномодифицированных сеток для решения задач сильноточной электроники // Вестник ЮУрГУ. 2012. № 40(299). С. 130 – 140. Свешников В.М., Рыбдылов Б.Д. О распараллеливании решения краевых задач на квазиструктурированных сетках // Вестник ЮУрГУ, серия ВМИ. 2013. Т.2, № 3. С. 63 – 72. Корнеев В.Д. Параллельное программирование в MPI // Новосибирск. ИВМиМГ СО РАН. 2002. Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. // Новосибирск: ИВМиМГ (ВЦ) СО РАН, 2001.
NVIDIA. CUDA C Best Practices Guide. 2012.
Studying the effectiveness of graphics accelerators paralleling solutions of three-dimensional boundary value problems on quasi-structured grids
I.A. Klimonov2, V.D. Korneev1, V.M. Sveshnikov1,2 When paralleling solutions of three-dimensional boundary value problems on quasi-structured grids by decomposition of the computational domain into subdomains, matched without the imposition of the most time consuming computational procedure is a solution of subproblems in subdomains. Using the parallelepiped quasi-structured grids makes possible to use for this purpose rapidly convergent method of alternating directions. The parallelization of the iterative process on subdomains is carried out in CPU and to solve the subproblems in this paper we propose to use the GPU graphics accelerators. Experimental studies of the use of graphics accelerators for solving subproblems by Peaceman – Rechford method is presented. We give experimental estimates of acceleration parallelization in hybrid computing environment CPU + GPU as compared to the calculations only on the CPU.
Keywords: boundary value problem, domain decomposition method, Poincare – Steklov equation, quasi-structured grids, Peaceman – Rechford method, graphic accelerators.