Параллельная реализация метода расщепления для системы из нескольких GPU с применением в задачах аэрогидродинамики Текст научной статьи по специальности «Математика»

Информационные технологии 246 Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 5 (2), с. 246-252

УДК 004.942

ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ МЕТОДА РАСЩЕПЛЕНИЯ ДЛЯ СИСТЕМЫ ИЗ НЕСКОЛЬКИХ GPU С ПРИМЕНЕНИЕМ В ЗАДАЧАХ АЭРОГИДРОДИНАМИКИ*

© 2012 г. С.Б. Березин1, И.С. Каргапольцев, Н.Д. Марковский2, Н.А. Сахарных1,2

1 Московский госуниверситет им. М.В. Ломоносова 2Компания NVIDIA Ltd, Московское представительство

nikolai.sakhamykh@gmail.com

Поступила в редакцию 10.09.2012

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

Ключевые слова: метод расщепления, численное моделирование течений, GPU.

Введение

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

Относительно недавно появились удобные инструменты для разработки приложений общего назначения для GPU - такие как CUDA, OpenCL, Parallel Nsight [1]. Эти средства значительно упростили разработку сложных приложений, эффективно использующих параллельную архитектуру GPU.

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

теров для получения результата достаточной точности за приемлемое время.

Решение уравнений Навье - Стокса на GPU

Течения жидкостей и газов в естественной среде (движение воздуха в земной атмосфере, воды в реках и морях, газа в атмосферах звезд и т.п.) и инженерных сооружениях (в трубах, каналах, струях, пограничных слоях около движущихся в жидкости или газе твердых тел, следах за такими телами и т.п.) могут быть описаны системой уравнений Навье - Стокса. Полная форма уравнений Навье - Стокса в трёхмерной декартовой системе координат для идеального газа с постоянной плотностью имеет вид

V • и = 0, (1)

Ян 1

— + u-Vu = -Vr + — V2u, (2)

dt Re

дТ 1 у -1

—— + u • V71 = —-— AT + ——- Ф, (3) dt Pr- Re y • Re

p = pRT, (4)

где u = (u. v. vi') - вектор скорости, p - давление,

p - плотность. Здесь уравнение неразрывности (1) выражает постоянство объема, (2) - закон сохранения движения, (3) - уравнение теплопроводности, и замыкает систему уравнение

Статья рекомендована к публикации Программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2012».

состояния (4). Число Рейнольса Яе = иЬ / V характеризует свойства течения: V = ц / р - кинематическая вязкость, ц - вязкость, а Ь и и -средняя скорость и размер области в конкретной задаче. В уравнениях также присутствуют число Прандтля Рг, газовая постоянная Я и удельная теплоемкость у. Ф - это диссипативная функция, характеризующая работу вязких сил:

Ф = Ф + Ф + Ф

х у г?

„, „I du У f ЭиУ f dw У du du dw du

Фх = 2 l — I +| — I +| — I +--+--,

dx J [dx J [ dx J dx dy dx dz

Ф = 2

y

\ 2

(du ^ ( дЛ

V

dy J I ay

du

2

( dw ^

2

/

du ]

2

Ф z = 21-I

z [dz J [dz

dy dw

+ du do + dw do dy dx dy dz'

\ +du c!w + do c!w dz J dz dx dz dy

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

Например, уравнение для х-компоненты (2)

ди ди ди ди --+ и--+ V--+ ^— =

д1 дх ду дг

дТ

1

■ +

dx Re

fd2u д2u д2uï

• + -

- + -

^dx2 dy2 dz2 расщепляется на три уравнения:

du du --+ u —

dt dx

dT 1 --+-

dx Re

f d2u ^

dx2

Подробный вывод уравнений, их описание и применение можно найти в [2].

Во многих случаях течение становится турбулентным - скорость, температура и давления начинают хаотично изменяться во времени и пространстве. В случае турбулентных течений численное решение системы (1)-(4) требует явного описания всего диапазона пространственных и временных масштабов, что возможно при числе узлов сетки пропорциональном Re94 [3]. Обычно турбулентность наступает при больших числах Re, порядка нескольких тысяч или десятков тысяч. Таким образом, практическое применение численного моделирования течений на основе системы уравнений Навье - Сто-кса сильно ограничено объемом памяти и производительностью существующих вычислительных систем. Существуют модели LES (large eddy simulation), основанные на воспроизведении наиболее значимых вихрей и параметризации оставшейся части спектра [3, 4], что позволяет снизить требования к разрешению сетки и вычислительным ресурсам, но требует выдвижения некоторых априорных утверждений о природе моделируемого течения, что не всегда представляеся возможными. В данной работе будет рассмотрен метод прямого численного решения системы уравнений Навье - Стокса без дополнительных упрощений.

Метод покоординатного расщепления

и соответствующий разностный метод

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

du du — + v— =

dt dy

du du --+ w —

dt dz

1

Re J_

" Re

fd 2u ^

dy2

fd 2u ^

dz2

(5)

(6)

(7)

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

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

Поскольку уравнения Навье - Стокса являются нелинейными, необходимо также проводить итерации для обновления нелинейных коэффициентов. В данной реализации проводятся итерации двух типов: глобальные и локальные. Глобальные итерации применяются для целого временного шага для всей системы, локальные - для каждого дробного шага, соответствующего одному из направлений. Общая схема алгоритма приведена на рис. 1. На каждом целом временном шаге система расщепляется по трем направлениям, и соответствующие системы уравнений последовательно решаются в глобальных итерациях. Число глобальных итераций выбирается так, чтобы численная ошибка, рассчитываемая как невязка в уравнении неразрывности, после каждого шага не превышала установленную величину. На каждом дробном шаге расщепления соответствующие уравнения решаются при использовании конечно-разностной схемы. На рис. 2 показана модель данных и расчетов для отдельного шага расщепления. Основное вычислительное ядро - это решатель наборов трехдиагональных систем. Полученные решения систем используются для обновления нелинейных коэффициентов. После этого проводится несколько локальных итераций для уменьшения численной ошибки.

+

2

+

п-1 шаг

п шаг

п+1 шаг

Рис. 1

Рис. 2

Подобная схема была использована для решения двумерных уравнений Навье - Стокса в работе [5]. В дальнейшем на основе этой идеи был реализован численный метод для трехмерного случая и динамическая система визуализации [6].

Реализация метода прогонок на одном GPU

В CUDA-реализации прогонки каждая нить решает ровно одну трёхдиагональную систему. За счет массивной параллельности нитей на

GPU и большого числа систем удается эффективно загрузить устройство. На рис. 3 показано распределение трехмерного массива на дробном шаге расщепления по отдельным нитям. Каждая нить обрабатывает отдельный столбец массива в определенном направлении. Код соответствующего ядра аналогичен CPU-версии (табл. 1).

Часто для достижения максимальной эффективности на GPU приходится менять структуру данных, а иногда и сам метод.

Анализ показал, что подобные прогонки вдоль оси Z работают медленнее, чем по X и Y.

Таблица 1

Реализация алгоритма прогонки на CUDA. Одна нить решает ровно одну систему

_device_void solve_tridiagonal(

FTYPE *a, FTYPE *b, FTYPE *c, FTYPE *d, FTYPE *x, int num, int id, int num_seg, int max_n )

{

get(c,num-1) = 0.0;

get(c,0) = get(c,0) / get(b,0); get(d,0) = get(d,0) / get(b,0); // прямой ход прогонки

for (int i = 1; i < num; i++) {

get(c,i) = get(c,i) / (get(b,i) - get(a,i) * get(c,i-1)); get(d,i) = (get(d,i) - get(d,i-1) * get(a,i)) /

(get(b,i) - get(a,i) * get(c,i-1));

}

get(x,num-1) = get(d,num-1); // обратный ход прогонки for (int i = num-2; i >= 0; i--) get(x,i) = get(d,i) - get(c,i) * get(x,i+1);

}

Это объясняется тем, что прогонки вдоль Ъ приводят к чтению каждой нитью последовательных значений в памяти. В результате, запросы нитей варпа к памяти не объединяются (ипсоа^се^ см. [1]).

акватории Белого моря (■«гЫ1е_8еа). В этом случае можно оценить насколько хорошо метод работает на реальных задачах со сложными границами и неравномерной загрузкой устройства. Тестирование производительности прово-

Рис. 3

Таблица 2

Влияние транспонирования на производительность (ср. время работы ядра в мс), NVIDIA Tesla C2050

Точность Одинарная Двойная

Ядро оригинал с трансп. ускорение Оригинал с трансп. Ускорение

ось X 523 528 1.0x 717 709 1.0x

ось Y 525 531 1.0х 694 686 1.0х

ось Z 1681 533 2.4х 1901 693 2.7х

трансп. 164 190

Всего 2729 1756 1.6х 3312 2278 1.5х

Таблица 3

Результаты тестов производительности на модели box_pipe, систем/с

Направ- CPU Tesla Tesla Tesla C2050

ление (4 ядра) C1060 C2050 vs CPU

ось X 1.55 5.61 17.73 11.4х

ось Y 2.17 5.36 18.11 8.3х

ось Z 2.34 5.64 18.12 7.8х

Все 1.96 5.30 16.09 8.2х

Таблица 4

Результаты тестов производительности на модели

white sea, систем/с

Направ- CPU Tesla Tesla Tesla C2050

ление (4 ядра) C1060 C2050 vs CPU

ось X 3.65 15.40 37.95 10.4x

ось Y 5.12 16.89 44.91 8.8x

ось Z 2.25 5.65 13.38 5.9x

Все 3.82 11.88 27.40 7.2x

Напротив, в прогонках вдоль X и У, соседние нити читают последовательные элементы в памяти, и все запросы объединяются. Оптимизировать прогонки по Ъ можно одним из двух способов: либо менять метод, либо менять шаблон доступа к данным. В данном случае реализован второй способ: входной массив данных транспонируется, а затем вместо Ъ-прогонки запускается У-прогонка. Эффективное транспонирование возможно за счет использования разделяемой памяти. В табл. 2 показано насколько быстрее работает улучшенный вариант с переупорядочиванием данных в одинарной и двойной точности.

Реализация метода расщепления тестировалась на двух задачах. В первой задаче моделировалось течение внутри прямоугольного канала со стенками, параллельными осям координат (Ьох_р1ре). Такой вид границ позволяет равномерно загрузить вычислительное устройство. Во второй задаче рассматривалась геометрия

дилось на одном GPU (Tesla C1060 или Tesla C2050), а также на одном многоядерном CPU (Intel Core i7, 4 ядра, 2.8GHz). На CPU метод был распараллен с помощью директив OpenMP и использовал все доступные ядра. В табл. 3 приведены результаты по каждому направлению расщепления и для всего метода в целом для теста box_pipe. В табл. 4 для сравнения приведены результаты тестов white_sea.

Реализация метода прогонок на нескольких GPU

Схема работы метода прогонок на нескольких GPU, подключенных к одному хосту, может состоять из следующих стадий:

- разделение данных между GPU;

- распараллеливание кода;

- синхронизация результатов между GPU.

Скорость обмена данными между GPU существенно ниже скорости глобальной памяти,

поэтому необходимо минимизировать количество пересылаемой информации. Разделение данных между GPU лучше всего провести вдоль направления X, поскольку все трехмерные массивы хранятся по строкам в виде одномерного вектора (i, j, k) ^ [dimY • dimZ • i + dimZ • j + k], и сетка естественным образом разбивается на блоки размера, кратного dimY • dimZ. В этом случае распараллеливание прогонок вдоль Y и Z тривиально и сводится к одновременному запуску ядер на нескольких GPU, работающих над соответствующими частями сетки (табл. 5).

Таблица 5

Одновременный запуск вычислительных ядер на нескольких GPU

for (int i = 0; i < nGPU; i++) {

cudaSetDevice(i); // Переключение устройства kernel<<<...>>>(devArray[i], ...); // Расчет своей

области сетки

}

Прогонка вдоль направления X разделяется между всеми GPU и выполняется последовательно, по цепочке: каждый GPU сначала ожидает данные, затем вычисляет часть результатов и отправляет их следующему GPU, используя cudaMemcpyPeer (рис. 4).

Рис. 4

Таким образом, с ростом числа GPU, прироста производительности вдоль оси X не будет, и, согласно закону Амдала, это послужит основным лимитирующим фактором произво-

дительности при дальнейшем масштабировании.

Тем не менее, другие части шага по X - вычисление производных и диссипативной функции - можно эффективно распределить между GPU, за счёт чего достигается хорошее масштабирование производительности на нескольких GPU, подключённых к одному хосту (рис. 5).

В силу того, что вычисление производных у границ требует доступа к граничным значениям частей сетки соседних GPU, необходима их эффективная синхронизация на каждом шаге метода. Синхронизацию такого рода можно реализовать с помощью функции прямого асинхронного копирования cudaMemcpyPeerAsync (или cudaMemcpyAsync при использовании UVA) (табл. 6). Также было проведено исследование влияния на производительность возможности прямого копирования данных между GPU в устройствах Tesla поколения Fermi, минуя оперативную память процессора.

1 □ 2 □ 4 □ 8 GPUs Рис. 5

Доступна реализация метода для нескольких GPU с использованием технологии MPI. Код можно использовать на нескольких нодах с различным количеством GPU на каждой из нод. Полный исходный код модели, использующей данную реализацию метода, можно найти на сайте проекта [7].

Таблица 6

Пример синхронизации между устройствами путем асинхронного копирования

// Массивы указателей с граничными полосами на каждом устройстве:

float **dev_halo_right; // правая синхронизируемая область

float **dev_halo_left; // левая синхронизируемая область

float **dev_data_right; // правая граничная область

float **dev_data_left; // левая граничная область

float **dev_data; // область, которую не требуется синхронизировать

cudaStream_t *devStream; // один поток на устройство

Таблица 6 (продолжение)

// Выделение и инициализация массивов и потоков на каждом устройстве. int haloSize = ...;

// Синхронизация правой области с левой на устройстве i+1. for( int i = 0; i < numDev - 1; i++)

cudaMemcpyPeerAsync(dev_halo_left[i+1], i+1, dev_data_right[i], i, sizeof(float) * haloSize, devStream[i]); // Синхронизация левой области с правой на устройстве i-1. for( int i = 1; i < numDev; i++)

cudaMemcpyPeerAsync(dev_halo_right[i-1], i-1, dev_data_left[i], i, sizeof(float) * haloSize, devStream[i]); // Ожидание окончания копирования:

for( int i = 0; i < numDev; i++) {

cudaSetDevice(i);

cudaStreamSynchronize(devStream[i]);

}

// Очередной полный пересчёт узлов в отдельных частях сетки.

for( int i = 0; i < numDev; i++) {

// Переключение индекса текущего GPU. cudaSetDevice(i);

kernel<<<grid, block, 0, devStream[i]>>>(dev_data[i], dev_data_left[i], dev_data_right[i], dev_halo_left[i], dev_halo_right[i]);

}

Заключение

В данной статье была представлена эффективная ОРИ-реализация метода покоординатного расщепления для решения полной системы уравнений Навье-Стокса в трехмерной области. Численный метод был протестирован на модельной задаче и продемонстрировал хорошие результаты по производительности.

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

Среди ближайших планов является импле-ментация параллельного алгоритма без ограничения масштабирования вдоль направления X. Для этого элементы сетки вдоль направления Ъ разбиваются на блоки ХУ, в которых выполнение прогонок вдоль X для одних блоков покрывается прогонками вдоль X и У для других. (Наличие локальных итераций вдоль каждого направления будет компенсироваться большим

количеством блоков). Ожидается, что эффективность данного алгоритма будет улучшаться с размером сетки.

Работа выполнена при финансовой поддержке РФФИ, проекты № 10-01-00288-а, № 09-07-00424-а.

Список литературы

1. NVIDIA Inc. NVIDIA CUDA programming guide, version 4.1. 2011.

2. Флетчер К. Вычислительные методы в динамике жидкостей. Том 2. М.: Мир, 1991.

3. Лапин Ю.В. Статистическая теория турбулентности // Научно-технические ведомости. 2004. 2(36). СПб.: Изд-во Политехнического университета.

4. Sagaut Pierre. Large eddy simulation for incompressible flows, 3rd edition. Springer, 2006.

5. Пасконов В. М., Березин С. Б. Неклассические решения классической задачи о течении вязкой несжимаемой жидкости в плоском канале // Прикладная математика и инфорика. №17: М.: МАКС Пресс, 2004.

6. Пасконов В.М., Березин С.Б., Корухова Е.С. Динамическая система визуализации для многопроцессорных систем с общей памятью и ее применение для численного моделирования турбулентных течений вязких жидкостей // Вестник Московского университета. Сер. 15. Вычислительная математика и кибернетика. 2007. C. 7-16.

7. Sakharnykh N., Berezin S., Paskonov V. Fluid solver based on Navier-Stokes equations using finite difference methods in areas with dynamic boundaries: URL: http://code.google.com/p/cmc-fluid-solver/

MULTI-GPU IMPLEMENTATION OF ADI METHOD WITH APPLICATIONS TO FLUID DYNAMICS S.B. Berezin, I.S. Kargapoltsev, N.D. Markovsky, N.A. Sakharnykh

Direct numerical simulation of turbulent flows requires highly detailed grids and huge computing resources. Modern GPUs have high floating-point instruction throughput and large memory bandwidth that makes them very suitable for computational fluid dynamics. Modeling complex flows demands extremely large grids. However, limitations on the device memory available on a single GPU makes it difficult to process complex models. Therefore, computations have to be distributed on a multi-GPU system. An efficient GPU-based parallel ADI algorithm and load balancing strategies for a multi-node system are proposed. A comprehensive performance analysis of the method is also presented for different input geometries.

Keywords: alternating direction implicit method, numerical simulation of flows, GPU.