CC BY
132
Уфа : УГАТУ, 2012
QjJ~AQ^lQj
___________________________________________Т. 16, №6(51). С. 108-112
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
УДК 519.3
Б. П. Рыбакин
РЕШЕНИЕ ТРЕХМЕРНЫХ ЗАДАЧ ГАЗОВОЙ ДИНАМИКИ НА ГРАФИЧЕСКИХ УСКОРИТЕЛЯХ
В работе представлены результаты расчетов трехмерных течений газа, полученных с помощью параллельных программ, разработанных для работы на многопроцессорных ЭВМ и графических ускорителях. Расчеты были проведены на CPU Intel Core I7 и GPU фирмы NVIDIA Tesla C1060 и GeForce 570 GTX. В алгоритме расчета использованы схемы высокой разрешающей способности типа TVD. Параллельные алгоритмы; газовая динамика; графические ускорители
Решение трехмерных задач газовой динамики представляет большой практический интерес. Эти решения используются при проектировании автомобилей и самолетов, исследованиях формирования галактик, взрывов сверхновых звезд, предсказаниях погоды, моделировании климата и т. д. В последние годы были разработаны эффективные разностные схемы высокого разрешения, которые позволяют проводить математическое моделирование указанных выше процессов. Но проведение трехмерных расчетов на подробных сетках требует больших вычислительных ресурсов. Поэтому появление новых технологий разработки программ для работы на графических ускорителях позволяет, в дополнение к существующим, использовать новые возможности для проведения интенсивных вычислений. Повышение производительности видеокарт привело к тому, что количество ядер в них достигло нескольких сотен. Например в Tesla C1060 240 ядер, а в новых графических ускорителях семейства Fermi (М2090) их количество достигает 512. На сегодняшний день существуют две основные технологии разработки программ для работы на GPU (Graphical Processor Unit) - CUDA и OpenCL. CUDA - это Compute Unified Device Architecture - технология программирования графических ускорителей, предложенная фирмой NVIDIA [1, 2]. Эта технология позволяет работать с языками программирования С и FORTRAN. OpenCL - Open Compute Language - совместная разработка ряда фирм, которая позволяет разрабатывать программы для графических ускорителей NVIDIA, ATI и обычных процессоров (CPU). Таким образом появляется возможность проведения интенсивных математических расчетов на графических ускорителях.
Модель программирования CUDA поддерживает четыре уровня абстракций: кооперативную организацию потоковых групп, общую память, барьерную синхронизацию внутри потоковых групп, а также координацию независимых потоков в grid (грид) [3]. GPU позволяют одновременно запускать на выполнение десятки тысяч потоков, причем минимальное количество потоков составляет один warp, который содержит 32 нити. Кроме того, графический ускоритель имеет свою собственную память, которая, в свою очередь, делится на глобальную, локальную, регистровую, память констант и текстур. Все это требует другого подхода к проектированию вычислительных алгоритмов и программирования.
ПОСТАНОВКА ЗАДАЧИ
Большой практический интерес представляет собой задача взаимодействия ударной волны с цилиндрической или сферической полостями, заполненными газом с другим значением плотности. Такие задачи характерны для моделирования сгорания топлива, испарения пузырьков и т. д. С другой стороны, такие задачи являются хорошими тестами для проверки различных алгоритмов, методов и программ. С практической точки зрения, такие задачи позволяют изучить процесс образования вихрей при контактном взаимодействии нескольких газов с различными плотностями и уравнениями состояния [4]. Систему уравнений газовой динамики, записанную в консервативной форме, можно представить в виде:
Э +f- <ри)=°-
Э ь, 0)
^ + — (pu, и, + РЪг1) = 0,
dt Эх, 11 1
Контактная информация: 8(495)939-40-06
^ + ^ [(е + Р)и, ] = 0. дt дх1
Уравнение состояния использовалось в виде: е = 1/2ру2 + е. В приведенных выше уравнениях (1) р - плотность, V^,{i = х, у, 2} - скорость, Р - давление, е - удельная внутренняя энергия, е - полная энергия.
Пусть область Я: {0 < х < 1, 0 < у < 1, 0 < 2 < < 1} представляет собой параллелепипед, заполненный покоящимся газом с плотностью р0 = = 1,0 и давлением р0 = 1,0. Внутри этого параллелепипеда находится сферическая область радиуса г = г0, с центром в точке (х0, у0, 20). Слева, в точке с координатами х = х*, располагается плоская ударная волна, которая движется влево. Начальные условия справа от ударной волны и вне сферического пузырька задаются следующим образом (р, р, и, V, w)т = (рй, рй, ий, ук, wh):. Внутри сферической полости давление и плотность равны р = рь =1,0; р = рь =0,1. За ударной волной величины определяются с помощью соотношений Ренкина-Гюгонио. Расчетный параллелепипед покрывался прямоугольной сеткой размерами от 64*64*64 до 1024x1024x1024.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ
Моделирование трехмерных газодинамических течений относится к числу сложных динамических процессов, которые предъявляют повышенные требования к используемым разностным схемам. Эти схемы должны максимально точно воспроизводить поведение вещества в окрестностях больших разрывов и достоверно описывать малые возмущения вдали от фронтов ударных волн. Для решения таких противоречивых требований используют разностные схемы, которые сочетают в себе повышенную разрешающую способность в областях малых возмущений и монотонность в областях сильных разрывов. К таким схемам принадлежат схемы типа ТУБ, БКО, WENO, РРМ. Такие нелинейные схемы второго порядка точности с ограничением общей вариации позволяют проводить расчеты ударных волн с большим разрешением и предотвращать нефизические осцилляции за их фронтами. Кроме того, такие алгоритмы должны хорошо распараллеливаться на многопроцессорные ЭВМ.
Для решения данной задачи была использована разностная схема повышенной разрешающей способности типа ТУБ. Условие общего ограничения осцилляций является нелинейным
условием устойчивости. В схеме ТУБ для дискретного решения его можно определить как
tv (ut)=X;
. , \uL - u‘\.
1+1 г+1 1
(2)
Условие (2) определяет меру общего количества осцилляций решения. Разностную схему записываем используя потоки второго порядка
точности К'+ш на границах ячеек [5]. Они получаются из потоков первого порядка точности
К
(i)t
1+1/2
superbee =
для противопотоковой схемы с помощью коррекции второго порядка точности. Поток первого порядка точности определяется из средних значений потока на границах ячеек. Коррекция второго порядка точности вводится для ограничения ложных осцилляций решения. Количество осцилляций на шаге по времени номер n + 1 не должно превосходить количество осцилляций на шаге по времени номер n.
TV(Un+1) < TV(Un).
Для минимизации осцилляций используются различные ограничители (limiters) - minmod, superbee, vanLeer. Ограничитель minmod позволяет выбрать наименьшее абсолютное значение между левой и правой коррекциями:
minmod(a, b) = [sign(a) + sign(b)]min(|a|, |b|).
Ограничитель superbee позволяет выбрать наилучшую коррекцию из:
minmod(a, 2b), if |a| > |b| minmod(2a, b), if |a| < |b|
Наиболее подходящим для данной схемы является ограничитель vanLeer из всех ограничителей, он позволяет находить гармоническое среднее значение между левыми и правыми коррекциями:
vanLeer(a, b) = 2ab .
a + b
По данной схеме была построена программа, которая использует алгоритм разделения по пространственным переменным (dimension splitting). Такая схема позволяет осуществлять эффективное распараллеливание как на многопроцессорных ЭВМ, так и на графических ускорителях. Решение, представленное в работе, получено с помощью схемы TVD с ограничителем vanLeer. Число Куранта выбиралось равным CFL = 0,5.
Для данного алгоритма было проведено тестирование [5, 6], которое показало, что предложенная схема и компьютерная программа позволяет с достаточной точностью решать газодинамические задачи, представляющие практи-
ческий интерес. Разностная схема и программа расчета была протестирована с помощью стандартных тестов задачи о распаде произвольного разрыва Сода и Лакса. Результаты численных расчетов показали хорошее совпадение с аналитическим решением [5].
лять на ОРИ, другая часть работала на центральном процессоре.
Рис. 1. График плотности на момент времени t = 5,021. Стрелками показаны векторы скорости для плотности р = 0,12
В данной работе представлены результаты расчетов взаимодействия ударных волн с областями газа пониженной плотности. Расчеты проводились на процессорах Intel Core I7 920 и графических ускорителях NVIDIA Tesla C1060 и GeForce 570 GTX. Графические ускорители Tesla C1060 и GeForce 570 GTX имеют разную архитектуру. В архитектуре Tesla есть 30 потоковых мультипроцессоров, каждый мультипроцессор имеет 8 ядер для работы с одинарной точностью, таким образом, общее количество ядер равно 240. Количество памяти, доступное для работы равно 4 Gb. Графический ускоритель NVIDIA GeForce 570 GTX имеет 480 ядер и 1.2 Gb памяти.
Версия программного обеспечения CUDA-Toolkit 4.0 дает возможность использовать сквозную адресацию памяти, которая физически принадлежит разным устройствам
Кроме того, можно запускать различные фрагменты программы на различные GPU. Программирование осуществлялось на языке Fortran фирмы PGI. Построение параллельного алгоритма для сжимаемых гидрогазодинамических течений было выполнено с использованием технологий OpenMP и CUDA для PGI Fortran. Для проведения расчетов на графических ускорителях программа расчетов была перестроена таким образом, чтобы часть подпрограмм направ-
Рис. 2. Распределение плотности на момент t = 5,G21 безразмерного времени
Рис. 3. Распределение энергии на момент времени t = 5,021
На графическом ускорителе проводились расчеты потоков на границах ячеек, вычислялась коррекция потоков и вычислялись ограничители. Архитектуры графических ускорителей Tesla и Fermi используют многопоточные процессоры SM, которые позволяют осуществлять распараллеливание с помощью технологии SIMT (Single Instruction Multiply Thread). Это позволяет создавать и управлять десятками тысяч конкурирующих потоков. Потоки можно организовывать в 1d, 2d и 3d потоковые блоки, каждый из которых работает на своем ядре. Все потоки выполняются независимо, используя потоковые регистры и локальную память. Для каждого потока и потокового блока вычислялись встроенные переменные: threadIdx
и blockIdx.
Рис. 4. Значение поля плотности на момент времени t = 12,51
На рис. 1 приведен график плотности на момент времени t = 5,G21. Внутри изоденсы p = = G,12 показано векторное поле для этой поверхности. На рис. 2 и З приведены сечения трехмерной расчетной области. На рис. 2 приведено распределение плотности, на рис. З -энергии. Разрешение фронта ударной волны равно размеру двух-трех пространственных ячеек. Таким образом, проведенное тестирование предложенных схем и программ трехмерного моделирования гидродинамических течений показывает, что данную схему можно применять для решения задач нестационарной газовой динамики, астрофизики и т. д.
Разрешение фронта ударной волны равно размеру двух-трех пространственных ячеек. Таким образом, проведенное тестирование предложенных схем и программ трехмерного моделирования гидродинамических течений показывает, что данную схему можно применять для решения задач нестационарной газовой динамики, астрофизики и т. д. [7, 8].
На рис. 4, 5 представлены графики плотности и энергии в моменты времени t = 12,51 и t = 15,G7 соответственно. На переднем фронте ударной волны, а также на некоторых внутренних поверхностях стрелками изображено векторное поле скорости. На рис. б приведен график плотности, энергии - рис. 7 на момент времени t = 2G,G57. Результаты проведенного трехмерного моделирования хорошо совпадают с результатами двумерных расчетов и позволяют проследить процесс деформирования пузырька под действием сильной плоской ударной волны.
Вычисления, проведенные на графическом ускорителе GeForce 57G GTX, дают ускорение по сравнению с одноядерным процессором на сетке 256*256x256 в 15,2 раза. При этом ос-
новная программа работает на центральном процессоре. На рис. 8 приведен график ускорения, полученный на 8 ядрах центрального процессора (СРи) для сеток 128*128*128 и 256*256*256.
Рис. 5. Распределение энергии на момент времени t = 12,51
Рис. б. Распределение плотности на момент времени t = 2G,G57
Рис. 7. Распределение энергии на момент времени t = 2G,G57
-
Асс Hers tion
wilh Up< пМЬ
-• -- --
;
у' / -12 3x12 8x1; 28
/.-■ -25 5x25 04 X со 56
1 2 3 4 5 6 7 8
Number of cores
Рис. 8. Ускорение вычислений на различных сетках
Анализ результатов, полученных при решении гидродинамических задач, практически невозможен без подробной графической визуализации полученных результатов. Это особенно верно для трехмерных расчетов. В данной работе графическая визуализация результатов компьютерного моделирования осуществлялась с помощью программы НОУК [9], любезно предоставленной В. Д. Горячевым, за что автор выражает ему глубокую признательность.
ВЫВОДЫ
В данной работе представлены разностная схема и программа для решения многомерных задач газовой динамики на графических ускорителях. Проведенное тестирование показало, что предложенные схемы и программы позволяют решать сложные газодинамические течения. Для этой цели использовалась схема ТУБ, обладающая высоким разрешением в областях ударных волн и зон больших градиентов. Вычисления были проведены на сетках размерами от 64^64x64 до 1024x1024x1024 узлов с помощью технологий ОрепМР и СЦОЛ. Ускорение вычислений для сетки размерами 128x128x128 и 512x512x512 приведено на рис. 8. Отметим, что используемый алгоритм достаточно эффективен для 8 ядерного процессора. Графические
ускорители помогают значительно ускорить вычисления.
СПИСОК ЛИТЕРАТУРЫ
1. NVIDIA CUDA. Programming Guide. 2011. http://developer. nvidia. com/ cuda-toolkit-40
2. CUDA FORTRAN. Programming Guide and Reference. 2011, http://www.pgroupe.com
3. Рыбакин Б. П. Параллельное программирование для графических ускорителей. М.: НИИСИ РАН. 2011, 261 с.
4. Kurganov A., Chi-Tien L. On the reduction of numerical dissipation in central-upwind schemes // Commun. Comput. Phys. February, 2007. Vol. 2, No. 1. P. 141-163.
5. Рыбакин Б. П., Шидер Н. И. Построение параллельных алгоритмов для решения задач гравитационной газовой динамики // Вычислительные методы и программирование. 2010. Т. 11. С. 388-394.
6. Рыбакин Б. П. Параллельная трехмерная TVD схема для решения задач гравитационной газовой динамики // Параллельные вычислительные технологии (ПаВТ'2009). Н.-Новгород, 2009. С. 673-679.
7. Rybakin B., Shider N. Computer modeling of multidimensional problems of gravitational gas dynamics on multiprocessor computers // Computer Science Journal of Moldova (CSJM). 2009. V. 17, № 1(49). P. 3-13.
8. Рыбакин Б. П. Решение III-d задач газовой
динамики на графических ускорителях. Сб. тр.
«Международная суперкомпьютерная конференция "Научный сервис в сети Интернет: экзафлопсное
будущее". 2011. C. 321-325
9. Горячев В. Д., Смирнов Е. М., Балашов М. Е. Визуализация результатов численного эксперимента при моделировании нестационарных течений с большим объемом данных в системе HDVIS // Научный сервис в сети Интернет: решение больших задач: Тр. Всероссийск. науч. конф. (22-27 сентября 2008 г., г. Новороссийск). М.: МГУ, 2008. C. 55-59.
ОБ АВТОРЕ
Рыбакин Борис Петрович, проф. каф. газовой и волновой динамики мех.-мат. факультета МГУ, вед. науч. сотр. НИИСИ РАН.
