CC BY
21
УДК 004.94
Компьютерное моделирование движения воздуха в условиях городской застройки с использованием NVIDIA CUDA1
А.Л. Архипов
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
Иваново, Российская Федерация E-mail: [email protected]
Авторское резюме
Состояние вопроса: Моделирование движения воздуха в условиях городской застройки является важной задачей, позволяющей определять направление распространения потоков воздуха в зависимости от расположения городских строений. Подобные задачи обычно решаются сеточными методами на многопроцессорных кластерах, позволяющих выполнять вычисления с большой скоростью. Однако стоимость такого компьютерного оборудования достаточно высока. Предлагается использовать в качестве вычислительного средства графический процессор. На сегодняшний день два основных игрока рынка графических процессоров - NVIDIA и AMD - предлагают технологии использования своих GPU для вычислений общего назначения (General-purpose computing on graphics processing units - GPGPU). Наиболее развитой технологией выполнения вычислений общего назначения на GPU является платформа NVIDIA CUDA.
Материалы и методы: Используется платформа NVIDIA CUDA, позволяющая выполнять вычисления на графических процессорах. Выбор этой платформы был обусловлен ее продуманным API и поддержкой самых распространенных на сегодняшний день видеокарт NVIDIA GeForce. Для использования графического процессора в качестве вычислительного средства адаптирована конечно-разностная схема вычислений.
Результаты: Предложен вариант решения задачи моделирования движения воздуха в условиях городской застройки с использованием графических процессоров современных видеокарт в качестве вычислительного устройства.
Выводы: Адаптированная конечно-разностная схема и алгоритм на ее основе позволили получить прирост скорости вычислений от использования графического процессора, по сравнению с вариантом алгоритма, ориентированным на работу с центральным процессором. В результате был получен значительный прирост скорости вычислений. Использование графических процессоров в качестве вычислительного средства и соответствующая адаптация алгоритмов позволяют снизить в 60-70 раз стоимость компьютерного оборудования, необходимого для вычислений, при решении задачи моделирования движения воздуха в условиях городской застройки.
Ключевые слова: численное моделирование, дифференциальные уравнения, уравнения Навье-Стокса, параллельные вычисления, NVIDIA CUDA.
Computer Simulation of Air Flows in Urban Areas Using NVIDIA CUDA
A.L. Arkhipov
Ivanovo State Power Engineering University, Ivanovo, Russian Federation E-mail: [email protected]
Abstract
Background: Modelling air movement in urban areas is an important problem, which allows determining the air flows depending on the location of city buildings. Such kind of problems is usually solved by mesh based methods on the multiprocessor supercomputer that allows performing calculations with great efficiency. However, this approach requires expensive computer hardware. In this article the usage of the GPU as a computational device instead of expensive multiprocessor computer is suggested. Nowadays, two main market players of graphic processors, NVIDIA CUDA and AMD, suggest the technologies of GPU usage for general purposes calculations. The more developing technology of general purposes calculations on GPU is NVIDIA CUDA platform.
Materials and methods: NVIDIA CUDA is used to calculate on the graphic processors. This choice was made on the basis of API and support of the most widely-spread NVIDIA GeForce video-cards. The finite-difference calculation scheme is adapted for using the graphic processors as a calculation means.
Results: In this article a way to solve a problem of the modeling the movement of air in urban areas is proposed. The GPU is used as a computational device instead of expensive multiprocessor computer.
Conclusions: The finite-difference scheme and the algorithm based on it allowed to get the velocity growth of calculation due to their implementation on the NVIDIA CUDA platform comparing to the algorithm which works with the central processor. As a result, the outstanding calculation velocity growth is received. The usage of graphic processors as a calculation means and algorithm adaptation allows to decrease the cost of computer equipment which is necessary for calculating, for solving the tasks of the movement modelling of air in urban areas up to 60 - 70 times.
Key words: numerical simulation, differential equations, Navier-Stokes equations, parallel computing, NVIDIA CUDA.
1 Работа выполнена при поддержке Министерства образования и науки РФ ГК№13 G25.31.0077.
Моделирование движения воздуха в условиях городской застройки может помочь найти ответы на вопросы, куда будут распространяться вредные выбросы того или иного промышленного объекта (например, электростанции или химического завода), где будет наибольшая концентрация вредного газа или дыма, исходящего из какого-либо источника техногенного или природного типа. Также такое моделирование может использоваться при решении вопроса о выборе места строительства экологически критичного промышленного объекта (например, АЭС или обычной тепловой электростанции).
Подобное моделирование обычно проводят с использованием параллельных вычислений для ускорения расчетов. В качестве вычислительных устройств применяются многопроцессорные компьютеры и вычислительные кластеры. Нами предлагается использовать в качестве вычислительных устройств графические процессоры (GPU), что, предположительно, позволит получить сравнимую с большими многопроцессорными компьютерами скорость вычислений при меньших материальных затратах.
На сегодняшний день два основных игрока рынка графических процессоров - NVIDIA и AMD - предлагают технологии использования своих GPU для вычислений общего назначения (General-purpose computing on graphics processing units (GPGPU)). Наиболее развитой технологией выполнения вычислений общего назначения на GPU является платформа NVIDIA CUDA [3].
В качестве дифференциальных уравнений, описывающих движение газа, в большинстве случаев используют уравнения Навье-Стокса:
Эр d(pu) d(pv) d(pw)
dt
dx
dy
dz
= 0;
d(Pu) , d(Pu +(Y- 1)Pem) , d(Puv)
dt
d(puw) dz
dx (d2(pu) dx 2
d2(pu) dy 2
dy
d2(pu) 1 dz2
= Qx
d(Pv) , d(Pvu) , d(Pv + (Y- 1)Pem) , d(Pvw)
dt
(
-v
dx
d2(pv) dx 2
d2(pv) dy 2
dy d2(pv) 1 dz2
dz
= Qy
d(pw) d(pwu) d(pwv) d(pw 2 + ( y- 1)pem)
dx
dy
dz
d2(pw)
d2(pw)
d2(pw) 1
ч dx dy dz
d(pE) + d(upE + u( y- 1)pem)
= Qz
dt
d(vpE -
dx
v (Y- 1)Pem)
5(w pE + w ( y- 1)pem)
dy dz
Qx = P9x + 2<»zPv; Qy = pgy - 2&zPu; Qz = P9z
= 0;
em = E--
w
2
R = (Y - 1)cv
где р - парциальная плотность; Т - температура; Е - удельная энергия; ет - внутренняя энергия 1 кг воздуха; u, v, w - скорости воздуха по Ox, Oy, Oz; R = 287 Дж/(кгК) - индивидуальная газовая постоянная для сухого воздуха; ®z - проекция вектора угловой скорости вращения Земли на вертикаль в данном месте; дХ, gy, gz - проекции ускорения свободного падения на оси Ox, Оу, От. gx = 0, gy = 0, gz = -9,81 м/с2; V - кинематическая вязкость воздуха; у = 1,4 - показатель адиабаты; оу = 717 Дж/(кгК) - теплоемкость воздуха при постоянном объеме; ор = 1006 Дж/(кгК) - теплоемкость воздуха при постоянном давлении.
Для удобства написания разностной схемы расчета следует представить дифференциальные уравнения Навье-Стокса в матричном виде.
dM
dt
M =
dFx dFv dFz
( d2V
dx dy dz
d2V
( p 1
pu
pv , Fx =
pw
,pe j
dx dy dz
pu 1
pu2 + (Y- 1)pem puv puw
u (pE + (Y- 1)Pem )
= Q,
fv =
( pv 1
pvu
pv2 + (y- 1)pem pvw
v (pE + (Y- 1)pem)
( 0 1 ( 0 pu pgx + 2<»zpv
pv , Q = pgy - 2rozpu pw pgz
.0 J I 0
E + (y- 1)em = E + RT = —
V =
Fz =
p w p wu pwv pw2 + (y- 1)pem w (pE +(Y- 1)pem)
w
■ + Yem =
w
cpT.
Если в качестве разностной схемы расчета выбрать обычную схему по методу конечных разностей, то придется использовать очень маленький шаг по времени, чтобы вычисления были устойчивыми. Лучше в качестве разностной схемы расчета выбрать такую схему, которая бы позволила использовать больший шаг по времени. Обычно это достигается за счет усложнения вычислений на каждом шаге интегрирования. Нами была выбрана сеточная схема, приведенная в [1]. Эта схема позволяет наиболее точно считать системы дифференциальных уравнений, содержащие в себе уравнения Навье-Стокса.
Схема вычислений использует дополнительные переменные и коэффициенты. т - шаг
2
по времени; п - момент времени; Л - шаг сетки по осям Ох, Оу, Oz; / - индекс по оси Ох, / = 2,...,М|-1, М1 - количество точек по оси Ох; 1 - индекс по оси Оу,} = 2,...,М2-1, N - количество точек по оси Оу; к - индекс по оси Oz, к = 2,...,М3-1, М3 - количество точек по оси Oz; 9 - сглаживающий коэффициент, 9 = 1,4.
Для получения значения матрицы М в следующий момент времени выполним вычисления.
МПм = мп1 к + т- М1 к,
где МП11к - состояние системы в точке сетки /,1,к в момент времени п.
Для получения коэффициента бМ1^ необходимо выполнить следующие вычисления [1].
Рхп 1 - Рхп 1 Руп 1 - Руп 1
/—,1 ,к '/+-,! ,к / 1—,к /,/'+—,к бмп1к = — 2' 21 + 1 2 1 2 +
2К
2h„
Fzn 1 - Fzn 1 i,j ,k— i,j ,k+—
----------2-------------г. і v • dVn, k і Qn, k;
2hz /, J ,k /, J ,k ’
Fxn 1,, = Fx
f
іі—, j,k 2
ML
/+—, j,k 2
MR
/+—, j ,k 2
іа 1
/і—, j,k 2
- MR1
ML 1
іі Туj ,k іі ^j,k
і 2 2 у
ML 1 = Mi, j,k і -/+г,J,^ 2
M. 1 = Mi і1, j,k /і—,J ,k
(Mx)i,j,k ;
2 ;
(Mx)iі1,j,k .
(Mx )i-1,j,k .
ml 1 = ml1 , = M/-1, j ,k 2
/—,J ,k /-1і—,j ,k 2
2^ 2
MR1 = MR1 1 = M-jk - ;
/—, J ,k /-1і—, j ,k 2
2 2
(Mx)/■ j,k = minmod| -M/ jM),
Miі1, j,k - Mi-1, j,k
minmod(ф1, ф2, ...) =
, 0(Mi, j,k - Mi-1,j,k));
min ф,, если ф, > 0V/,
І
max ф,, если ф, < 0V/,
f
І і—, j ,k 2
= тах
Г
дFv
0, в противном случае;
f f
ML
І і—, j ,k 2
дМ
,r
дFv
MR1
І і—, j ,k і 2,J’ у
дМ
где г - радиус распространения возмущений.
Для расчета подобных аэродинамических моделей максимальная скорость распространения возмущений а в точке пространства определяется следующим образом.
І і—, j,k 2
= тах
= тах
i, j і—,k 2
І і—, j ,k 2
І, j і—,k 2
л
L R R
і C 1 , uR 1 і C 1
І і х, J,k i і x, J,k І і ~, j ,k
2 2 2
\
, „L R , „R
і C 1 , v 1 і C 1
i, J і o,k /, J і ~,k І, j і ~,k
2 2 2 )
І, j ,k і— 2
= тах
w
І, j ,k і— 2
л
L R R
і CL 1, wR 1 і CR 1
І, j ,k і x І, j ,k і x І, j ,k і x
2 2 2 )
где с - скорость звука в данной точке пространства.
0, 1,к =лМУ- 1)ет .
Для расчета влияния вязкости можно
применить явный метод конечных разностей.
1/п 2^п + 1/п //'-1,1',к - 2//, 1,к + //+1,1 ,к
dVinj ,k =■
h2
Vn - 2Vn і Vn Vn - 2Vn і Vn V,j-1,k 2/i,j,k і Vi,j!1,k , Vi,j,k-1 2/j,j,k і Vi,j,kі1
h?
h2
Значения скоростей и, V, w, температура и давление в точке вычисляются по формулам.
pu pv
u =—, v = —, p p
e = pE = (pu )2
em = =
p
f
w =
pw
(pv)2 і (pw)2
2p2
T =
P = (Y-1)
pE -
(pu)2 і (pv)2 і (pw) 2p
После получения значений всех переменных на текущем шаге по времени можно переходить к вычислению следующего шага. Условие устойчивости схемы.
1
-ma,x (yfij
/,J,k V ’
vi, j,k і wi, j,k
■ Ci,j,k I <
8
тіп(Ьх ,Ьу ,Ь2)
Реализация. Нами предлагается использовать вычисления на графическом процессоре (СРи) для реализации алгоритма расчета вышеописанной разностной схемы.
СРи представляет собой потоковый процессор, состоящий из большого числа однотипных ядер. Количество ядер достигает нескольких сотен, что много больше, чем количество ядер СРи (центрального процессора). Использование такой мультиядерной архитектуры позволяет достичь существенного возрастания скорости вычислений. Однако у СРи есть и ограничения: потоковый процессор обеспечивает только параллелизм данных, но не параллелизм задач. То есть можно выполнять одну и ту же операцию над большим количеством данных одновременно, но нельзя запустить одновременно несколько разных задач. При разработке алгоритмов для выполнения на СРи нужно учитывать эту особенность.
Рассмотренная выше вычислительная схема позволяет создать алгоритм с использованием параллелизма данных: в каждой точке сетки выполняются одинаковые вычисления для нахождения значений переменных в сле-
v
2
дующий момент времени. Таким образом, данная вычислительная схема может быть использована для выполнения расчетов на СРЫ.
Синхронизация вычислений для каждого момента времени была организована следующим образом. после вычисления значений всех переменных на СРЫ управление передавалось СРи, результаты вычисления обрабатывались, после чего происходил переход к следующей итерации алгоритма. Без такой синхронизации вычисления потеряли бы устойчивость. Понятно, что такая синхронизация, как и вообще сама организация вычислений на СРЫ, требует определенных издержек - часть вычислительного времени тратится на передачу данных от СРЫ к СРи и обратно, а также на обслуживание процессов, связанных с самими вычислениями на СРЫ. Поэтому прирост скорости вычислений относительно СРЫ заметен тем больше, чем большего размера сетка обрабатывается.
Каждый объект, находящийся в расчетной области, препятствует движению воздуха. Для учета этого на его гранях задаются соответствующие граничные условия. В качестве таких граничных условий был применен следующий прием. внутри объектов скорость воздуха обнуляется.
Вычисления расчетной схемы на СЫРА были организованы таким образом, что каждое ядро СРЫ выполняло расчет состояния одной ячейки расчетной области. В него передавались переменные р, ри, ру, р^, рЕ. Для уменьшения задержек, связанных с обращением к глобальной памяти видеокарты, данные ячеек, необходимых для выполнения вычисления, кэшировались в регистрах. Шаг по времени т, показатель адиабаты у и кинематическая вязкость воздуха V были заданы константами. На каждом шаге сначала вычислялись новые значения переменных, а затем накладывались граничные условия. в ячейках, примыкающих к границам расчетной области, и ячейках, принадлежащих объектам-препятствиям, значения переменных соответствующим образом корректировались. На данном этапе, чтобы излишне не усложнять алгоритм, границы объектов задавались непосредственно в коде программы.
Исходный код программы был составлен таким образом, что его можно было скомпилировать как для расчетов с использованием СЫРА, так и для расчетов с использованием СРЫ. В варианте расчетов на СРЫ вычисления производились в цикле по всем ячейкам в один поток. Такая организация исходного кода позволила выполнить сравнение времени выполнения этих двух вариантов программы и оценить полученное ускорение.
Выводы и результаты. Были проведены вычисления 100000 шагов по времени на сетке разной размерности. Шаг сетки составлял 1 м, шаг по времени - 0,0001 с. Время работы программы в секундах занесено в таблицу.
Размер сетки GPU CPU
50x60x60 180 10000
64x64x64 344 20900
100x100x100 1331 70900
Описанная вычислительная схема в варианте для GPU обсчитывалась в среднем в 60 раз быстрее однопотокового варианта на CPU. Соответственно, чтобы получить такую же скорость вычислений на многопроцессорном кластере, он должен содержать не менее 60 процессорных ядер. Стоимость компьютера с видеокартой NVIDIA GeForce будет в 60-70 раз меньше стоимости такого многопроцессорного кластера. Однако производимые вычисления потребовали для себя существенного объема памяти видеокарты, поэтому размер расчетной области был ограничен - максимальный размер куба составил 100x100x100 ячеек. В расчетной области были заданы два препятствия в виде параллелепипедов: один высотой б0 ячеек (метров), другой - 32 ячейки (метра). Внутри них скорость воздуха обнулялась. На левой грани куба расчетной области была задана начальная скорость входящего потока воздуха по осям (x, y, z) - (б, б, 0) м/с. Ось x направлена слева направо, ось y - сверху вниз. Программа позволяет просматривать срезы расчетной области по оси z (рис. 1, 2).
Рис. 1. Скорость воздуха по осям x, y на уровне z = 20
Рис. 2. Скорость воздуха по осям x, y на уровне z = 33
По осям х и у максимальная величина скорости составила 5,1, по оси z - 1,2.
Дальнейшее развитие предполагает введение дополнительной субстанции в состав воздуха (например, дыма), картину распространения которой надо будет определить. Также планируется задавать форму и расположение объектов в расчетной области динамически, а не непосредственно в коде программы, как это было сделано в данном слу-
чае. Также планируется изменить отображение результата вычислений: вместо горизонтального среза расчетной области показывать трехмерный вид расчетной области. Кроме того, планируется увеличить размерность сетки, реорганизовав код таким образом, чтобы сократить требуемый для вычислений объем памяти.
Список литературы
1. Kurganov A., Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations // J. Comput. Phys. - 160 (2000). - Р. 241-282.
2. Kurganov A., Levy D. A third-order semi-discrete central scheme for conservation laws and convection-diffusion equations // SIAM J. Sci. Comput. - 22 (2000). - Р. 1461-1488.
3. NVIDIA CUDA — неграфические вычисления на графических процессорах [Электронный ресурс]. - Режим доступа: http://www.ixbt.com/video3/cuda-1.shtml]
4. Википедия: CUDA [Электронный ресурс]. - Режим доступа: http://ru.wikipedia.org/wiki/CUDA]
References
1. Kurganov, A., Tadmor, E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. J. Comput. Phys., 160 (2000), рр. 241-282.
2. Kurganov, A., Levy, D. A third-order semi-discrete central scheme for conservation laws and convection-diffusion equations // SIAM J. Sci. Comput., 22 (2000), рр. 1461-1488.
3. NVIDIA CUDA - negraficheskie vychisleniya na graficheskikh protsessorakh [General-purpose computing on graphics processing units]. Available at: http://www.ixbt.com/video3/cuda-1.shtml]
4. Vikipediya: CUDA [Electronic Resource]. Available at: http://en.wikipedia.org/wiki/CUDA]
Архипов Ален Леонидович,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», программист ИВЦ, аспирант кафедры высокопроизводительных вычислительных систем, e-mail: [email protected]
б
