CC BY
123
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Том 152, кн. 1 Физико-математические пауки 2010
УДК 519.684:532.546
РЕШЕНИЕ ЗАДАЧ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ С ПРИМЕНЕНИЕМ ТЕХНОЛОГИИ
NVIDIA CUDA
Д.Е. Демидов, А.Г. Егоров, А.Н. Нуриев
Аннотация
На примере модельной задачи о развитии неустойчивости Кельвина Гельмгольца рассмотрены преимущества технологии NVIDIA CUDA при решении задач гидродинамики численными методами. Рассмотрены как сеточные методы, так и бессеточный. Приводится краткое описание технологии CUDA.
Ключевые слова: вычислительная гидродинамика, параллельные вычисления, технология CUDA.
Введение
В настоящее время графические процессоры (GPU) являются оптимальной по соотношению цена производительность параллельной архитектурой с общей памятью. Действительно, видеокарта NVIDIA GTX 285. доступная на рынке по цене около 15 ООО руб.. имеет пиковую производительность 950 GFLOPS. тогда как один из последних четырехядерных центральных процессоров (CPU) Intel Core2 Extreme QX9775 всего 102 GFLOPS и продается по цене около 48 ООО руб. Причина такой производительности видеокарт заключается в том. что они изначально были сконструированы для одновременного применения одной и той же шейдерной функции к большому числу пикселей, или. другими словами, для высокопроизводительных параллельных вычислений.
До недавнего времени использовать вычислительные мощности видеокарт было крайне неудобно, так как программисту приходилось выражать свою задачу в графических терминах. Ситуация изменилась в 2007 г.. когда компания NVIDIA выпустила технологию CUDA программно-аппаратную архитектуру для вычислений на графических процессорах. Несмотря на свою молодость. CUDA уже приобрела широкую популярность в научных кругах. Вычислительные задачи, реализованные на CUDA. получают ускорение в десятки и сотни раз в таких областях, как молекулярная динамика [1, 2]. астрофизика [3. 4]. медицинская диагностика [5. 6] и др.
В настоящей работе рассматривается возможность применения технологии CUDA к задачам вычислительной гидродинамики. В качестве тестовой выбрана классическая двумерная задача о развитии неустойчивости в свободном сдвиговом слое между противоположно направленными потоками вязкой жидкости. Изначально бесконечно тонкий вихревой слой на границе потоков вследствие неустойчивости Кельвина Гельмгольца распадается на систему вихревых сгустков (рис. 1). С течением времени ширина вихревого слоя растет за счет спаривания вихревых сгустков [7. 8]. Качественный статистический анализ данного процесса требует выполнения большого объема вычислительных экспериментов, что оказывается возможным с использованием технологии CUDA.
% 4* * ^ X ^ ^ 9 Л - * • * * # \ * * .
(a) t = 1
^ j/ • . ' ** *« • •
(b) t = 3
/*- * • . • (%J
(с) t = 5
Рис. 1. Развитие неустойчивости Кельвипа-Гельмгольца в свободном сдвиговом слое
Для решения соответствующей задачи вычислительной гидродинамики широко используются как сеточные, так и бессеточные методы. Наиболее популярные из них рассмотрены в настоящей работе. Это. с одной стороны, геометрический и алгебраический многосеточные методы, а также метод редукции с использованием быстрого преобразования Фурье при решении задачи на сетке. Бессеточные методы представлены простейшим вариантом метода точечных вихрей. Во всех случаях при решении задачи на GPU получено значительное ускорение (десятки и сотни раз) по сравнению с CPU.
1. Программная модель CUDA
CUDA это среда разработки приложений, которая предоставляет в распоряжение программиста вычислительные мощности современных графических карт NVIDIA. В основе программной модели CUDA лежат три ключевые абстракции: иерархия групп потоков, разделяемая память и барьерная синхронизация, реализуемые в виде минимального набора расширений языка С/С---.
Приложения CUDA используют как центральный, так и графический процессоры. Фактически видеокарта является высокопроизводительным сопроцессором. Параллельные участки кода выполняются на видеокарте в виде вычислительных ядер. Ядро это функция, которая вызывается из основной программы на CPU
threadldx: 01234567
пппп
float х = input[threadldx]; float у = func(x); output[threadldx] = y;
Рис. 2. Ядро выполняется массивом потоков, каждый из которых имеет уникальный идентификатор
и выполняется на видеокарте. Для выполнения каждого ядра запускается большое число одновременных потоков. Укажем два ключевых, на наш взгляд, различия между программными потоками на центральном процессоре и на видеокарте. Во-первых. десятки тысяч потоков CUDA могут быть созданы всего за несколько процессорных циклов. Благодаря этому расходы на управление потоками значительно ниже, чем на CPU. и даже легковесные ядра, состоящие всего из нескольких строчек кода, могут быть чрезвычайно эффективными. Во-вторых, переключение между контекстами потоков происходит мгновенно н используется для маскировки задержек при обращении к памяти. При этом если один из потоков простаивает в ожидании операции чтения/записи или синхронизации, он заменяется другим потоком, выполняющим арифметические операции. Это означает, что для эффективной загрузки видеокарты число потоков должно многократно превышать число процессоров.
Итак, каждое ядро CUDA выполняется массивом параллельных потоков. Все потоки выполняют один и тот же программный КОД. II для того чтобы они могли обрабатывать свой участок данных, каждому из них назначается уникальный идентификатор. Этот идентификатор используется для вычисления адресов памяти и организации ветвлений. Например, па рис. 2 представлено ядро, выполняемое восемью потоками. Сначала каждый из них использует свой идентификатор для считывания элемента из входного массива. Затем каждый поток передает считанное значение в функцию и записывает результат в выходной массив.
Это пример так называемой неограниченно параллельной (eriibarassirigly parallel) задачи, когда каждый поток может выполнять свою подзадачу независимо от соседей. В реальных задачах требуется организация межпоточного взаимодействия: потоки должны, во-первых, совместно использовать результаты, чтобы избежать лишних вычислений: во-вторых, разделять операции доступа к памяти, чтобы снизить требования программы к пропускной способности канала данных. В программной модели CUDA взаимодействие допускается только между потоками внутри небольших групп, или блоков. Каждый блок это одно-, двух- или трехмерный массив потоков фиксированного размера. При запуске ядра потоки организуются в одно- или двухмерную решетку (grid) блоков. Потоки внутри блока могут взаимодействовать посредством разделяемой памяти и выполнять барьерную синхронизацию. Организация потоков в блоки позволяет программам прозрачно масштабироваться при переходе на видеокарту с большим числом потоковых процессоров.
Поясним последнее утверждение. Благодаря тому, что блоки потоков независимы. они могут выполняться в любом порядке на любом нз доступных процессоров. Допустим, ядро состоит из 8 блоков потоков. На видеокарте с двумя мультипроцессорами каждый мультипроцессор должен будет выполнить по 4 блока. А па видеокарте с 4 мультипроцессорами каждый из них должен будет выполнить
Глобальная память, Локальная память
Мультипроцессор
Блок
двойной
Текстурный
процессор
Разделяемая
память,
Регистры
Потоковые
процессоры
Рис. 3. Архитектура GPU 10-й серии
Различные типы памяти, доступные программам па CUDA
Табл. 1
Вид Область видимости Скорость доступа Объем Использование
Регистровая Поток Высокая 16384 регистра па МП Локальные переменные потока
Локальная Поток Низкая Ограничен объемом глобальной памяти Локальные переменные. пе уметившиеся в регистрах
Разделяемая Блок Высокая 16 КВ Организация межпоточпого взаимодействия
Глобальная Программа Низкая До 4 ГВ Хранение данных, обмен с CPU
по 2 блока, то ость задача будет выполнена в два раза быстрее. Такое масштабирование происходит автоматически без каких-либо изменений в программе. Поэтому однажды написанная программа на CUDA может выполняться на самых различных по производительности устройствах от ноутбуков до высокопроизводительных серверов.
На рис. 3 представлена архитектура видеокарт NVIDIA 10-й серии (например. GTX 280). Видеокарта содержит 240 потоковых процессоров, каждый из которых может одновременно выполнять до 96 потоков. Потоковые процессоры сгруппированы в 30 .мультипроцессоров (МП), каждый из которых содержит 8 потоковых процессоров, блок двойной точности, текстурный процессор и блок регистровой и разделяемой памяти.
Потоки CUDA имеют доступ к нескольким областям памяти, которые различаются областью видимости, скоростью доступа и объемом (табл. 1). Регистровая память используется для хранения локальных переменных потока. Переменные, которые не уместились в регистрах, хранятся в так называемой локальной памяти, которая, несмотря на название, находится вне чипа и потому имеет низкую скорость доступа. Разделяемая память используется для организации взаимодействия между потоками в блоке. Все потоки ядра имеют доступ чтения/записи к глобальной памяти, которая располагается в памяти видеокарты. Область видимости глобальной памяти все приложение: ее содержимое не изменяется между запус-
ками различных ядер. Кроме того, центральный процессор также имеет доступ к глобальной памяти, поэтому глобальная память используется для обмена данными между СРи и СРи.
Рассмотрим модель исполнения кода в С1ГОА. Каждый поток выполняется потоковым процессором: блок потоков выполняется одним мультипроцессором: каждый мультипроцессор может выполнять несколько блоков. Число блоков, которое может выполняться на одном мультипроцессоре, ограничено разделяемыми ресурсами. то есть объемом регистровой и разделяемой памяти.
2. Сеточные методы
Перейдем к решению тестовой задачи о развитии неустойчивости Кельвина Гельмгольца. При использовании подхода Эйлера течение вязкой жидкости в свободном сдвиговом слое описывается уравнениями Навье Стокса. В терминах функции тока (ф) — завихренности (ш) они имеют вид:
Аф = ш, (1а)
дш . (дф дф \
Областью течения является полоса —Н < у < Н. На границах у = ±Н, моделирующих бесконечность, ставятся условия ф = ш = 0. В продольном направлении ф и ш считаются периодическими функциями с периодом Ь. Начальное условие соответствует бесконечно тонкому сдвиговому слою:
Ь = 0: ш = 2(1 + С(х)Жу).
Здесь £ - функция Дирака, С _ малое случайное возмущение.
Дискретизация уравнения (1) проводится на равномерной прямоугольной сетке методом конечных разностей со вторым порядком точности по пространственной переменной. Для решения дискретной задачи используется полуиоявиая схема с поэтапным решением уравнений для завихренности и функции тока, предложенная в [9]. При этом уравнение переноса завихренности (1Ь) решается явно методом Рун-ге Кутты четвертого порядка, так что для вычисления значения завихренности на новом временном слое требуется 4 шага. После каждого шага необходимо обращать оператор Лапласа для определения функции тока. Схему решения можно
представить следующим образом:
ш __ шП
1 /9 +(ц"-Уь)ц"=£А„ьЛ ф1 = А^ъ (2а)
т/2
ш __шп
2 /9---Ь (111 • V;, VI = ьАНШ1, ф2 = Д^о^2, (2Ь)
т/ 2
шз__шп
—^----1- (и2 • УЛ.)^2 = еДл.^2, фз = А^0шз, (2с)
к4 = —т [(из •Уь) шз — еА^шз], (2(1)
<^п+1 = д (—шп + + 2и>2 + шз) + — фп+1 = ДЙ д(^>п+1. (2е)
Верхними индексами здесь обозначены значения переменных на соответствующих временных итерациях; нижними индексами обозначены промежуточные значения. Нижний индекс Н обозначает оператор, дискретизированный на сетке с шагом Н. Запись А-0(-) обозначает решение задачи Пуассона с нулевым граничным условием Дирихле для заданной правой части.
Табл. 2
Видеокарты, использовавшиеся в вычислительных экспериментах
Модель Стоимость (руб-) Число МП Пиковая производительность (GFLOPS) Пропускная способность памяти (ГБ/с)
8500 GT 1500 2 43 13
8600 GTS 3000 4 139 32
GTX 260 7000 27 804 112
Использование метода четвертого порядка позволяет обеспечить устойчивость данной схемы при разумном выборе шага по времени т ~ h. Наиболее дорогостоящей операцией, очевидно, является решение уравнения Пуассона для функции тока. Здесь необходимы алгоритмы, обладающие высокой степенью параллельности. например, редукция с использованием быстрого преобразования Фурье [10] нлн более общий многосеточный метод [11] в его геометрическом либо алгебраическом варианте. Производительность этих методов оценивалась с использованием трех видеокарт NVIDIA 8500 GT, 8600 GTS и GTX 260 (табл. 2).
2.1. Геометрический многосеточный метод. Многосеточные методы (ММ) считаются оптимальными для решения эллиптических дифференциальных уравнений [11] и основываются на использовании сеточной иерархии и операторов перехода от одной сетки к другой. Основная идея заключается в ускорении сходимости итерационного метода, лежащего в основе ММ, за счет коррекции решения, полученного на точной сетке, с помощью решения уравнения на грубой сетке.
Геометрический многосеточный метод (ГММ) основывается на предопределенной сеточной иерархии. При решении сеточного уравнения Au = / на каждом уровне иерархии задается сетка П;, оператор задачи A;, оператор продолжения P; : u;+i ^ u; и оператор сужения R; : u; ^ u;+1.
На каждом уровне l иерархии один шаг ГММ выполняет следующие действия:
• текущее решение сглаживается несколькими шагами релаксации;
• вычисляется невязка r; = / — A; u;, которая затем сужается в правую часть более грубого уровня: /;+1 = R;r;;
• если уровень l + 1 - последний в иерархии, уравнение A;+1u;+1 = /;+1 решается прямым методом; иначе на уровне l + 1 рекурсивно выполняется один шаг многосеточного метода с нулевым начальным приближением u;+1 = 0;
• решение на уровне l корректируется: u; ^ u; + P;ui+1;
•
Наиболее дорогостоящей операцией здесь является релаксация. Она занимает до 60% от общего времени счета, поэтому реализация этой процедуры особенно важна. В настоящей работе используется красно-черная релаксация Гаусса Зойдоля, так как она естественным образом распараллеливается.
2.2. Алгебраический многосеточный метод. Алгебраический многосо-точиый метод (АММ) [11] является одним из наиболее эффективных методов решения разреженных неструктурированных систем линейных уравнений большой размерности. В отличие от геометрического многосеточного метода (ГММ),
1 /* Умножение разреженной матрицы, Л на вектор х. */
2_global____void spmv(int m , int к, float *Adat, int *Aidx, float *x, fl oat *y) {
3 /* Каждый поток определяет номер своей строки матрицы, */
4 int row = blockldx.x * block І )ім і. х + throadldx.x;
<
6 int j, col;
7 float val;
8 float buf = 0.0;
В /* Цикл no ненулевым, элементам, строки */
<
11 col = Aidx[row + j * m];
12 val = Adat[row + j * m];
13 /* Проверка, что ненулевые элементы в строке не кончились * /
<
15 buf += val * x[col];
16 }
17 v[row] = buf;
18 }
IB }
Рис. 4. Умножение разреженной матрицы в формате ELL па вектор
АММ но тробуот постановки задачи на сетке, а работает непосредственно с разре-
Au = /
Алгоритм АММ состоит из двух основных этапов, примерно равных по числу операций: этап настройки метода и этап решения. Этап настройки заключается в конструировании проблемно-зависимой иерархии и включает в себя выбор вложенных подмножеств П;+1 С П; исходных перемеиных По, построение операторов продолжения Р; : П;+1 ^ П; и операторов та грубых сетках A;+1 = P;TAiPi. Эта часть алгоритма АММ основана только на алгебраической информации отно-A
Au = /
ких правых частей, то этап настройки достаточно выполнить один раз. Алгоритм этапа решения АММ практически совпадает с алгоритмом ГММ.
Нами АММ применяется для решения задачи Пуассона в (2). Этап настройки в классическом варианте АММ практически невозможно распараллелить, поэтому он полностью выполняется на CPU. Однако это не сказывается на скорости решения задачи (1), так как этап настройки достаточно выполнить один раз.
Этап решения был полностью реализован на GPU. Основным алгоритмическим элементом здесь является умножение разреженной матрицы на вектор. К этой операции сводятся практически все компоненты метода. Операции умножения разреженной матрицы на вектор с использованием CUDA посвящено несколько статей, например [12, 13]. Более всего па производительность влияет выбранный формат хранения разреженной матрицы. Одним из наиболее эффективных форматов является ELL [12], обеспечивающий когерентный доступ к памяти (coalesced memory access) соседними потоками [14]. В этом формате разреженная матрица размерности M х N с максимальным числом K ненулевых элементов в строке хранится в плотном массиве Adat размера M х K, где строки с числом ненулевых эле-K
хранятся в массиве Aidx размерности M х K. Оба массива размещаются в памяти по столбцам. На рис. 4 приведен пример ядра, выполняющего умножение разреженной матрицы на вектор. Каждый поток в ядре выполняет умножение одной строки матрицы на вектор и заполняет соответствующий элемент выходного массива.
Табл. 3
Решение задачи с применением сеточных методов
Устройство Время решения Ускорение
Геометрический мпогосеточпый метод
CPU AMD At lilonfi і 3200 9 120 с
GPU 8500 GT 1 824 с 5х
GPU 8600 GTS 608 с 15х
GPU GTX 260 190 с 48 х
Алгебраический мпогосеточпый метод
CPU AMD Лі lilonfi і 3200 14 009 с
GPU 8500 GT 2 838 с 5х
GPU 8600 GTS 1 044 с 14х
GPU GTX 260 311 с 45 х
Метод редукции
CPU AMD Лі lilonfi I 3200 4 560 с
GPU 8500 GT 798 с 6х
GPU 8600 GTS 264 с 17х
GPU GTX 260 87 с 52х
Другими часто выполняемыми операциями являются вычисление скалярного произведения векторов, определение максимального элемента в массиве, суммирование элементов массива. Библиотека с открытым исходным кодом CUDPP [15] позволяет эффективно выполнить эти операции на видеокарте.
2.3. Метод редукции. Метод редукции с применением быстрого преобразования Фурье [10] хорошо известен н применяется для нахождения решения простейших сеточных эллиптических уравнений в прямоугольнике. Основными алгоритмическими компонентами метода являются быстрое преобразование Фурье и обращение трехдиагональных матриц. Нами использовалось быстрое преобразование Фурье из библиотеки CUFFT [16]. входящей в среду разработки CUDA. Библиотека позволяет выполнить прямое и обратное построчное преобразование Фурье вещественной матрицы. При обращении системы трехдиагональных матриц каждый поток решал одну систему линейных уравнений вида
u1 — un — ° ui-1 (2 + 4sin ki)ui + ui+1 — h
2.4. Оценка эффективности. В табл. 3 приведены результаты решения задачи (1) с применением схемы (2) описанными выше методами. Задача решалась на сетке 1024 х 1025 в области (—п, п) х (—п, п), было выполнено 4000 шагов по времени. Как видно, достигаемые ускорения для каждого из графических процессоров практически не зависят от выбранного метода. Причем даже наименее производительная видеокарта 8500 GT позволяет достичь по крайней мере пятикратного ускорения. Для данной задачи оптимальным оказался метод редукции. Однако он является узкоспециализированным методом, в то время как многосеточные методы применимы к гораздо более широкому классу задач.
Как видно из табл. 3. ускорение, получаемое сеточными методами, прямо пропорционально пропускной способности видеокарт. Это же подтверждает и рис. 5.
3. Метод частиц
Рассматриваемая задача может быть сформулирована и на основе подхода Лагранжа [17]. При этом поле завихренности представляется в виде конечного набора
50
40
о
■| 30
и
"и
о
о
<
20
10
60 80 100 120 0 200 400 600 800
GB/sec GFLOPS
(а) (Ь)
Рис. 5. Зависимость ускорения сеточных методов (а) и метода частиц (Ь) от пропускной способности и производительности видеокарт соответственно
точечных вихрей (частиц). Развитие поля завихренности во времени определяется движением этих частиц. В пренебрежении диффузией закон движения вихрей
определяется системой обыкновенных дифференциальных уравнений
^Xj(#) = ^K(xj(#);xj(#)), 0 < i < N, (За)
2п?
х*(0) = ^еь (3b)
где Xj(t) - позиция i -й частицы в момент времени t, а функция K задает попарное гидродинамическое взаимодействие между двумя частицами. С учетом периодичности эта функция имеет следующий вид:
К(х; у) =---- ^^ d = х - у.
cosh d2 — cos di у — sin di J
Для учета диффузии в таком описании задачи может использоваться метод случайного блуждания (random walk) [18]. При этом на каждом временном шаге позиции частиц изменяются на случайную величину, распределенную по нормальному закону с дисперсией а2 = 2ет, где е - коэффициент диффузии, ат - величина шага по времени. Для генерации последовательности псевдослучайных чисел был использован алгоритм Морсонна Твистора (Morsonne Twister) [19]. Так как этот алгоритм носит чисто последовательный характер, каждый поток использовал собственную копию генератора.
Заметим, что задача (3) - это типичная задача о взаимодействии N тел. Такие задачи встречаются в астрономии, молекулярной динамике, гидродинамике и идеально подходят для реализации на GPU [20].
В табл. 4 приведены результаты решения задачи (3) методом частиц для числа вихрей N = 1024, N = 4096 и N = 8192. Как видно, карты 8500 GT и 8600 GTS сразу же выходят на максимальную производительность. Карта GTX 260 в первых двух экспериментах используется не полностью, с чем н связано относительно низкое ускорение. Так, в первом эксперименте при размере блока в 256 потоков
Табл. 4
Время решения задачи (3). В каждом из экспериментов использовался шаг по времени т = 10-
Устройство Время решения Ускорение
N = 1024, Ттах = 20
CPU AMD At hlonfi і 3200 5 353 с
GPU 8500 GT 118 с 45 х
GPU 8600 GTS 36 с 149х
GPU GTX 260 30 с 179х
N = 4096, Tmax = 5
CPU AMD At hlonfi 1 3200 21 156 с
GPU 8500 GT 408 с 51х
GPU 8600 GTS 129 с 164х
GPU GTX 260 37 с 565х
N = 8192, Tmax = 1
CPU AMD At hlonfi 1 3200 16 025 с
GPU 8500 GT 313 с 51х
GPU 8600 GTS 105 с 152х
GPU GTX 260 20 с 783 х
задействованы всего 4 мультипроцессора из 27 доступных. Поэтому неудивительно. что в этом эксперименте ускорения для карт 8600 GTS (4 мультипроцессора) и GTX 260 практически совпадают.
На порядок большие ускорения, достигаемые методом частиц по сравнению с сеточными методами, объясняются тем. что скорость вычислений здесь лимитируется производительностью карты, а не ее пропускной способностью (рис. 5).
4. Статистическая обработка результатов
Полный анализ полученных результатов выходит за рамки настоящей работы. Приведем лишь график зависимости ширины сдвигового слоя l от времени (рис. 6). В каждом вычислительном эксперименте величина l определялась как
l{t) = у - J и{х, у,1)у2с1хс1у.
Как видно на рис. 6. отдельные вычислительные эксперименты демонстрируют нерегулярное поведение, особенно при турбулизации сдвигового слоя при t > 1.
l
Для ее нахождения было проведено несколько сотен вычислительных экспоримон-
l
линиой. В начальный период времени соответствующая зависимость хорошо описывается известной формулой ламинарного режима I = aJ-net (пунктирная линия). При переходе к турбулентному режиму эта зависимость становится линейной. Заметим, что для проведения такого количества вычислительных экспериментов без применения технологии CUDA вместо двух дней потребовалось бы несколько месяцев.
5. Обсуждение
Как видно, любой из рассмотренных нами методов получает значительный прирост в производительности при реализации на GPU. Ускорение сеточных методов
Рис. 6. Зависимость ширины I сдвигового моя от времени осредненная по 750 реализациям (жирная сплошная линия). Тонкими линиями представлены зависимости Ь(Ь) для двух отдельных экспериментов. Области I. II и III соответствуют ламинарному, переходному и турбулентному режимам
прямо пропорционально пропускной способности памяти видеокарты. Это объясняется тем. что отношение числа арифметических операций к числу обращений к глобальной памяти в этих методах невелико:
Таким образом, сеточные методы не используют полную вычислительную мощность видеокарты, так как процессоры в основном простаивают в ожидании очередной порции данных. В методе частиц отношение P ^ 1, благодаря чему ускорение при реализации на GPU прямо пропорционально пиковой производительности видеокарты.
Пожалуй, самым серьезным недостатком расчетов на GPU является то. что вычисления с одинарной и двойной точностью существенно различаются по производительности. Это объясняется тем. что в современных видеокартах на каждые восемь потоковых процессоров приходится один модуль двойной точности. Однако этот недостаток оказывается несущественным для сеточных методов, ограниченных пропускной способностью памяти. Опыт показывает, что при использовании двойной точности производительность сеточных методов снижается всего в 1.3
1.5 раз.
Работа выполнена при финансовой поддержке РФФИ (проект X- 08-01-00548а).
Summary
D.E. Demidov, A.G. Egorov, A.N. Nuriev. Application of NVIDIA CUDA Technology for Numerical Solution of Hydrodinamic Problems.
Advantages of NVIDIA CUDA technology are presented on the example of classical problem of evolution of Calvin Helmliolz instability. Finit.e-difference as well as mesliless methods are considered. Short observation of CUDA technology is given.
Key words: computational fluid dynamics, parallel computations, CUDA technology.
Литература
1. Stone J.E.E., Phillips J.C.C., Freddolino P.L.L. et al. Accelerating molecular modeling applications with graphics processors // J. Comput. Cliem. 2007. V. 28, No 16.
P. 2618 2640.
2. Van Meel J.A., Arnold A,, Frenkel D. et al. Harvesting graphics power for MD simulations j j Mol. Simulat. 2008. V. 34. No 3. P. 259 266.
3. Zwart S.F.P. , Belleman R.G., Geldof P.M.. High-performance direct gravitational N-body simulations on graphics processing units // New Astronomy. 2007. V. 12. No 8.
P. 641 650.
4. Harris C., Haines K., Staveley-Smith L. GPU accelerated radio astronomy signal convolution // Exp. Astron. 2008. V. 22. No 1. P. 129 141.
5. Muyan-Ozcelik P., Owens J.D., Xia J., Samant S.S. Fast deformable registration on the GPU: A CUDA implementation of demons // Computational Science and its Applications: Intern. Conf. 2008. P. 223 233.
6. Noiil P.B., Walezak A., Hoffmann K.R. et al. Clinical Evaluation of GPU-Based Cone Beam Computed Tomography // Proc. of High-Performance Medical Image Computing and Computer-Aided Intervention (HP-MICCAI). 2008.
7. Winant C.D., Browand F.K. Vortex pairing: the mechanism of turbulent mixing-layer growth at moderate raynolds number // J. Fluid Mecli. 1974. V. 63, No 2. P. 237 255.
8. Aref H., Siggia E.D. Vortex dynamics of the two-dimensional turbulent shear layer // J. Fluid Mech. 1980. V. 100, No 4. P. 705 737.
9. Weinan E., Liu J.-G. Vorticity boundary condition and related issues for finite difference schemes // J. Comput. Pliys. 1996. V. 124. P. 368 382.
10. Hockney R. W. A fast direct solution of Poisson’s equation using Fourier analysis // J. ACM. 1965. V. 12, No 1. P. 95 113.
11. Trottenberg U., Oosterlee C., Schillier A. Multigrid. London: Acad. Press, 2001. 631 p.
12. Bell N. Efficient sparse mat.rix-vect.or multiplication on CUDA: NVIDIA Technical Report NVR-2008-004 / N. Bell, M. Garland: NVIDIA Corporation, 2008.
13. Baskaran M.M. Optimizing sparse mat.rix-vect.or multiplication on GPUs: IBM Research Report. RC24704 (W0812-047) / М. M. Baskaran, R. Bordawekar: IBM, 2009.
14. NVIDIA CUDA Programming guide. NVIDIA Corporation, 2009. Version 2.2.
15. Harris M. CUDA data parallel primitives library. URL: http://gpgpu.org/
developer/cudpp.
16. CUDA CUFFT Library. NVIDIA Corporation, 2009. Version 2.2.
17. XoKim P., Иствуд Д. Численное моделирование методом частиц. М.: Мир, 1987. 639 с.
18. Marsden J.E., Chorin A.J. A Mathematical Introduction to Fluid Mechanics (Texts
in Applied Mathematics). Springer, 1993. 169 p.
19. Matsumoto М., Nishimura T. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator // ACM Trans. Model. Comput.. Simul. 1998. V. 8, No 1. P. 3 30.
20. Nyland L., Harris М., Prins J. Fast. N-body simulation with CUDA // GPU Gems 3 / Ed. by H. Nguyen. Addison Wesley Professional, 2007. P. 677 695.
Поступила в редакцию 28.12.09
Демидов Денис Евгеньевич кандидат физико-математических паук, старший научный сотрудник Казанского филиала Межведомственного суперкомпыотерпого центра РАН.
Е-шаП: ddemiduvQksu.ru
Егоров Андрей Геннадьевич доктор физико-математических паук, заведующий кафедрой аэрогидромехапики Казанского государственного университета.
Е-шаП: Andrey.EguruvQksu.ru
Нуриев Артем Наилевич бакалавр механики, магистрант кафедры аэрогидромехапики Казанского государственного университета.
Е-шаП: аНет.501 Qlist.ru
