CC BY
41
УДК 519.688 Дата подачи статьи: 03.10.16
DOI: 10.15827/0236-235X.114.143-148 2016. Т. 29. №° 4. С. 143-148
ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ КАК СРЕДСТВО ПОВЫШЕНИЯ ЭФФЕКТИВНОСТИ РЕШЕНИЯ ЗАДАЧ ВЫЧИСЛИТЕЛЬНОЙ
АЭРОДИНАМИКИ
С.А. Буренков, аспирант, burenkovsaJ@mPei■ru; О.Ю. Шамаева, к.т.н., доцент, shamayevaoy@mpei.ru (Национальный исследовательский университет «Московский энергетический институт», ул. Красноказарменная, 14, 111250, г. Москва, Россия)
Предлагаются программные средства повышения эффективности решения класса аэродинамических задач на высокопроизводительных вычислительных системах кластерного типа. Программные средства реализованы на основе параллельно-последовательного алгоритма решения разреженных систем линейных алгебраических уравнений большой размерности методом обобщенных минимальных невязок.
В статье дается теоретическая оценка трудоемкости разработанного алгоритма, а также выполнено сравнение результатов теоретической оценки и результатов вычислительного эксперимента.
Программные средства реализованы для модели вычислений с распределенной памятью и использованы для решения модельной задачи обтекания профиля воздушным потоком. Приведены результаты вычислительного эксперимента по расчету плотности и давления обдувающего профиль воздушного потока на высокопроизводительном кластере, иллюстрирующие эффект от применения параллельной программы. В статье исследованы характеристики ускорения параллельного решения модельной задачи, позволившие определить порог размерности решаемых систем уравнений, ниже которого организация параллельных вычислений отрицательно сказывается на эффективности использования вычислительных ресурсов. Исследовано также влияние степени разреженности матриц коэффициентов на ускорение параллельной реализации.
Область применения результатов настоящей работы - не только вычислительная гидро- и газодинамика, но и энергетика, биология, химия, фармакология, астрофизика и др.
Ключевые слова: вычислительная аэродинамика, задача обтекания, параллельные вычисления, модель передачи сообщений MPI.
Развитие современного авиастроения требует решения разнообразных сложных задач, таких как обеспечение безопасности авиаперевозок, автоматизация систем управления полетами, создание двигателей, работающих на альтернативном топливе, и многих других. Помимо этого, не перестает расти потребность в увеличении грузоподъемности летательных аппаратов, скорости их полета и тяги силовых установок. Решение большинства из перечисленных проблем можно получить, прямо или косвенно прибегая к помощи аэродинамических механизмов расчета разного рода течений.
Аэродинамические характеристики обдувающего какой-либо объект потока можно определить в ходе натурного эксперимента, который проводится в специальной аэродинамической трубе. Эта процедура достаточно трудоемкая, дорогостоящая, а ее точность существенно зависит от условий проведения. Тем не менее до конца прошлого века аэродинамические испытания проводились исключительно средствами натурного эксперимента.
В настоящее время решение задачи аэродинамического обтекания можно получить с использованием методов математического моделирования. Такой подход позволяет с большой степенью точности определить все необходимые характеристики без использования громоздких стендов и проведения дорогостоящих натурных испытаний. Однако моделирование аэродинамических процессов относится к классу сложных задач [1, 2] в силу
высокой трудоемкости, больших объемов необходимой памяти, а также наличия специфичных для аэродинамики свойств. Например, расчет аэродинамических характеристик турбулентного течения с числом Рейнольдса порядка 105 в канале методом прямого численного моделирования на архитектуре производительностью 1 ТФлоп/с займет 122 года [1]. Поэтому процессы моделирования требуют использования специальных численных методов, параллельных технологий и программ, ориентированных на высокопроизводительную вычислительную технику, такую как суперкомпьютеры и кластерные системы.
На сегодняшний день известны несколько подходов к моделированию аэродинамических течений, среди которых метод прямого численного моделирования (Direct Numerical Simulation, DNS), метод моделирования крупных вихрей (Large Eddy Simulation, LES), метод моделирования разделенных вихрей (Detached Eddy Simulation, DES) и др. [1, 2]. Многие модели, такие как LES и DES, основываются на существовании средней скорости потока и среднего отклонения от нее. Поэтому реализация осредненных моделей проигрывает по точности моделирования модели DNS, но требует значительно меньшей вычислительной работы. Высокая вычислительная трудоемкость препятствует использованию DNS при решении комплексных задач аэродинамики, таких как обтекание планера самолета воздушным потоком, даже с использованием высокопроизводительной суперкомпьютер-
ной техники. На основании оценки развития вычислительных ресурсов, выполненной в работе [1], применительно к решению задачи обтекания среднего гражданского самолета или автомобиля расчет методом DNS можно будет производить за приемлемое время к концу XXI века. Более того, данная оценка является неоправданно оптимистичной, так как предполагает рост производительности аппаратных средств в два раза каждые пять лет. Появившиеся не так давно подходы к моделированию турбулентных течений (LES, DES и др.), как правило, менее трудоемки, но их использование увеличивает погрешность получаемого решения. К тому же каждый из них имеет свои преимущества, недостатки и ограниченную область применения [2]. Тем не менее использование DNS оправдано при моделировании обтекания объектов несложной геометрии.
Постановка задачи аэродинамического обтекания
В работе рассматривается задача аэродинамического обтекания профиля воздушным потоком. Предлагаются способы повышения эффективности ее решения на высокопроизводительной вычислительной технике с использованием механизма параллельных вычислений. Требуется определить основные аэродинамические характеристики - скорость, давление, число Маха и др. при обдувании воздушным потоком профиля RAE2822, изображенного на рисунке 1.
Рис. 1. Профиль RAE2822 Fig. 1. RAE2822 profile
Определить указанные характеристики потока можно, решив систему дифференциальных уравнений Навье-Стокса:
dW 1
р-= Я.-grad(р) + цAW + grad(divW) .(1)
А 3
Получить аналитическое решение системы (1) для большинства практических задач не представляется возможным, поэтому систему Навье-Стокса решают численными методами, для чего вокруг обдуваемого профиля строится расчетная область, изображенная на рисунке 2, которая разбивается неравномерной прямоугольной расчетной сеткой, как показано на рисунке 3.
После дискретизации системы (1) дальнейшее решение задачи сводится к системе линейных алгебраических уравнений (СЛАУ) большой размерности, матрица которой обладает свойствами разреженности и нерегулярности портрета, отсутствием диагонального преобладания и симметрии.
Рис. 2. Расчетная область для профиля RAE2822
Fig. 2. Computational domain for RAE2822 profile
Рис. 3. Фрагмент расчетной сетки профиля RAE2822 Fig. 3. A detail of RAE2822 profile computational grid
Результаты решения задачи обтекания могут быть представлены графически в виде распределения основных аэродинамических характеристик вокруг обдуваемого профиля. На рисунке 4 изображены распределения давления и скорости в воздушном потоке.
Организация параллельно-последовательных вычислений
На основе анализа методов решения СЛАУ с учетом перечисленных свойств матрицы коэффициентов [3, 4] предпочтение было отдано проекционному методу обобщенных минимальных невязок (General Minimal Residual, GMRES) [5], на основе которого в работе предложена и реализована параллельная модификация GMRES для кластерных архитектур.
Метод обобщенных минимальных невязок относится к классу методов проекционного типа и основывается на проектировании в подпространствах Крылова [5-7].
На рисунке 5 изображена структурная схема параллельного алгоритма решения СЛАУ по методу GMRES, вычисления в которой представлены в виде совокупности взаимодействующих последовательно выполняемых процессов. Для организации взаимодействия используется коммуникационная библиотека передачи сообщений (Message Passing Interface, MPI) [8, 9].
б)
Рис. 4. Распределение основных аэродинамических характеристик обдувающего профиль RAE2822 воздушного потока: а) давления, б) скорости
Fig. 4. Distribution of the main aerodynamic characteristics of airflow blowing RAE2800 profile: а) pressure, б) velocity
Главный процесс рассылает подчиненным заданные размерность подпространства Крылова m, требуемую точность решения eps, максимальное количество итераций maxNumbIter, а также соответствующие им части матрицы коэффициентов A и вектора правой части Ь (блоки 2 и 3). В качестве начального приближения к решению выбирается нулевой вектор. Производится предварительная оценка невязки начального приближения г, и вычисляется ее норма betta: i-й процесс производит вычисление фрагмента вектора невязки г(,), после чего фрагменты собираются главным процессом (блок 5).
Далее происходит итерационный процесс вычислений (блоки 6-11), завершающийся при достижении заданной точности или при выполнении максимального числа итераций. На каждой итерации процессы совместно формируют базисные вектора w, которые представляются в матричном виде H (блок 8). Полученная матрица H приводится к форме Хессенберга (блок 9) и используется в качестве матрицы коэффициентов треугольной СЛАУ H ■ у = g (блок 10). Решение полученной треугольной СЛАУ используется при уточнении решения x исходной системы уравнений, полученного на предыдущей итерации. Наконец, по вышеописанному принципу вычисляются новая невязка системы г и ее норма betta (блок 11). Отметим,
что наиболее вычислительно трудоемкие фрагменты алгоритма выполняются параллельно (блоки 5, 8 и 11).
Предложенный алгоритм обладает свойством масштабируемости и предполагает равномерную загрузку всех доступных вычислительных узлов кластера как исходными данными, так и непосредственно вычислительной работой [10].
Оценка трудоемкости параллельного алгоритма
Оценим вычислительную сложность реализованного параллельного алгоритма на множестве вычислительных узлов p при следующих предположениях. Пусть п - размерность решаемой СЛАУ, m - размерность подпространства Крылова, L - количество выполненных итераций, - длительность аддитивной арифметической операции, Р - длительность мультипликативной арифметической операции, ё - доля ненулевых элементов в матрице коэффициентов, которая определяется как
ё = ——— , где N - количество всех коэффициен-
N
тов матрицы, а N0 - количество коэффициентов, равных нулю. Тогда трудоемкость Т представленного в работе алгоритма может быть оценена с помощью формулы
Т (п, р, ё, т, Ь, Б, Р ) =
= O
( 3 2 Д
— L ■ d-\m + il-S + —L ■ d-\m + il-P p p
(2)
В работе выполнено сравнение теоретической оценки трудоемкости алгоритма с экспериментальными данными. На рисунке 6 представлены графики зависимости времени решения СЛАУ с помощью разработанной модификации GMRES от количества используемых вычислительных узлов. Для этой же системы уравнений был проведен вычислительный эксперимент, результаты которого отмечены на графике пунктиром.
Отметим, что время решения СЛАУ в ходе эксперимента немного больше предсказанного оценкой независимо от количества используемых вычислительных узлов. Этот эффект объясняется тем, что оценка (2) не учитывает временные затраты на чтение исходных данных и затраты, необходимые для осуществления коммутационных обменов между исполнителями.
Результаты вычислительного эксперимента
Вычислительный эксперимент проводился на вычислительном кластере из 32 процессоров AMD Opteron 6272 с тактовой частотой 2.1 ГГц, соединенных сетью Gigabit Ethernet, под управлением
1. Начало
2. Чтение входных данных: размерность подпространств (т), точность (eps), максимальное количество итераций (maxNumЫter), матрица коэффициентов вектор правой части (Ь)
1
3. Распределение фрагментов матрицы A и вектора b по исполнителям
1
4. p = m numbIter = 0
5. Вычисление невязки и ее нормы
r(0)=b(0)-A(0)x(0)
r(1)=b(1)-A(1)x(1)
r(n-1 )=b(n-1 )-A(n-1 )x(n-1 )
betta(0)=||r(0)||2
betta(1)=||r(1)||2
betta(n-1)=||r(n-1)||2
X
r, betta=max(betta(i))
12. Вывод полученного решения и числа итераций
7. numbIter = numbIter + 1 G = betta * e1
13. Конец
8. Построение базиса в подпространстве Крылова, упаковка базисных векторов в матрицу H
w(0) w(1) w(n-1 )
9. Приведение матрицы H к форме Хессенберга
10. Решение треугольной СЛАУ: H*y = g
11. Получение точного решения, вычисление невязки этого решения и ее нормы
x(0)=x(0)+sum(y(0)*v(0))
ж
x(1)=x(1)+sum(y(1)*v(1 ))
1
x(n-1)=x(n-1 )+sum(y(n-1)*v(n-1 ))
r(0)=b(0)-A(0)x(0)
r(1)=b(1)-A(1)x(1)
betta(0)=||r(0)||2
"С
ж.
r(n-1 )=b(n-1 )-A(n-1 )x(n-1 )
betta(1)=||r(1)||2
X
betta(n-1)=||r(n-1)||2
и
betta=max(betta(i))
7
Рис. 5. Обобщенная схема параллельно-последовательного алгоритма GMRES Fig. 5. The generalized scheme of GMRES parallel-sequential algorithm
ОС Linux Kernel 3.12.13. Для оценки эффективности параллельного решения использовалось понятие ускорения параллельного решения, определяемого как S(p)=Ti/Tp, где T - время последовательного решения; Tp - время параллельного решения с использованием p однородных вычислительных
узлов. На рисунках 7 и 8 приведены зависимости достигаемых ускорений решений от количества вычислительных узлов. При анализе зависимостей (рис. 7) можно отметить, что для разреженной системы уравнений небольшой размерности (п = 103) ускорение с ростом количества используемых вы-
H
16000 14000 12000 10000 8000 6000 4000 2000 0
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39
-Теоретическая оценка
- Эксперимент
Рис. 6. Сравнение теоретической оценки времени работы разработанной модификации GMRES с экспериментальными данными
Fig. 6. Comparison of theoretical estimates of the time of the developed GMRES modification and the experimental data
8 12 16 20 24 28 32 Количество вычислительных узлов
Рис. 7. Зависимости ускорений решения СЛАУ от количества вычислительных узлов для разных размерностей систем
Fig. 7. Dependences of SLAE decision accelerations on the number of computing nodes for different system dimensions
20 18 16 14 е12 ни10
& 8 о
о 6
^ 4
2 0
8 12 16 20 24 28 Количество вычислительных узлов
- d=0.01
d=0.001
d=0.0001
Рис. 8. Зависимости ускорений решения СЛАУ от количества вычислительных узлов для разной степени разреженности матриц коэффициентов
Fig. 8. Dependences of SLAE decision accelerations
on the number of computing nodes for systems with different degree of sparsity of coefficient matrixes
числительных узлов убывает. Это явление закономерно и объясняется тем, что выигрыш от приме-
нения распределенных вычислении не компенсирует затраты на их организацию. Для систем уравнении большей размерности (n > 104) ускорение решения растет с ростом числа вычислительных узлов. Таким образом, применение параллельных технологий не всегда оправдано; в частности, для рассмотренного примера целесообразно использование последовательного варианта метода GMRES при малых размерностях СЛАУ, а уже для систем размерности n > 104 - параллельного варианта.
На временные характеристики решения СЛАУ, а также на эффект от применения параллельных вычислений оказывает влияние и степень разреженности d матрицы коэффициентов. На рисунке 8 изображена зависимость ускорения решения СЛАУ от количества используемых вычислительных узлов кластера для разной степени разреженности матрицы коэффициентов размерностью n = 105. Отметим, что степень разреженности матрицы оказывает существенное влияние на ускорение процессов решения задачи на кластерных системах, начиная с числа используемых вычислительных узлов p > 20: чем более плотная матрица подлежит обработке, тем выше ускорение решения.
Таким образом, чем больше размерность системы и сложнее структура матрицы ее коэффициентов, тем больший эффект можно получить от применения параллельных технологий.
В заключение отметим, что использование высокопроизводительной вычислительной техники без масштабируемых параллельных программ, к сожалению, не приводит к росту эффективности процессов решения сложных задач, а в ряде случаев даже увеличивает время решения. Созданное параллельное программное приложение, использующее механизм параллельных вычислений, позволяет существенно уменьшить временные затраты на проведение вычислительных экспериментов и решать сложные аэродинамические задачи.
Литература
1. Беляев К.В., Гарбарук А.В., Стрелец М.Х., Шур М.Л., Спаларт П.Р. Опыт прямого численного моделирования турбулентности на суперкомпьютерах // Суперкомпьютерные дни в России: сб. тр. конф. молод. ученых. М., 2016. С. 357-364. URL: http://russianscdays.org/files/pdf16/357.pdf (дата обращения: 02.10.16).
2. Гарбарук А.В., Стрелец М.Х., Шур М.Л. Моделирование турбулентности в расчетах сложных течений. СПб: Изд-во Политехнич. ун-та, 2012. 88 с.
3. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. М.: Мир, 1991. 370 с.
4. Писсанецки С. Технология разреженных матриц. М.: Мир, 1988. 398 с.
5. Крукиер Л.А., Пичугина О.А., Чикина Л.Г. Обзор методов подпространства Крылова // Современные проблемы математического моделирования: сб. тр. XIV Междунар. конф. Ро-стов-н-Д: Изд-во ЮФУ, 2011. 400 с.
6. Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. Новосибирск: Изд-во РГГУ, 2000. 70 с.
7. Жуков В.Т., Новикова Н.Д., Феодоритова О.Б. Сдвиговая стратегия в обобщенном методе минимальных невязок. М.:
2
0
4
36
"n=100000"
4
32
36
Изд-во ИПМ РАН, 2009. 28 с.
8. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб: БХВ-Петербург, 2002. 608 с.
9. Гергель В.П. Теория и практика параллельных вычислений. М.: БИНОМ. Лаборатория знаний, 2007. 423 с.
10. Буренков С.А., Шамаева О.Ю. Библиотека эффективных функций для решения класса аэродинамических задач на кластерных архитектурах // Актуальные вопросы развития инновационной деятельности в новом тысячелетии: сб. тр. XIV Междунар. науч.-практич. конф. 2015. N° 3 (14). С. 8-11.
Software & Systems Received 03.10.16
DOI: 10.15827/0236-235X.114.143-148 2016, vol. 29, no. 4, pp. 143-148
PARALLEL CALCULATIONS AS A WAY TO INCREASE THE EFFICIENCY OF COMPUTING AERODYNAMIC PROBLEM SOLUTIONS
S.A. Burenkov1, Postgraduate Student, burenkovsa@mpei.ru
O.Yu. Shamayeva1, Ph.D. (Engineering), Associate Professor, shamayevaoy@mpei.ru 1 National Research University "MPEI", Krasnokazarmennaya St. 14, Moscow, 111250, Russian Federation
Abstract. The article proposes software to increase the efficiency of solving aerodynamic problems in highperformance cluster systems. The software is implemented on the basis of the proposed parallel-sequential algorithm for solving sparse systems of linear algebraic equations of high dimensionality using the method of generalized minimal residual.
The article provides a theoretical estimate of the developed algorithm complexity, as well as a comparison of the results of theoretical estimates and computational experiment results.
The software is implemented for a distributed memory computation model and is used to solve a model problem of an aerodynamic flow of a profile with an air flow. There are the results of a computational experiment on calculating the density and pressure of the blowing profile air flow in the high-performance cluster, which illustrate the effect of using a parallel program. The article analyzes the characteristics of model problem parallel solution acceleration, which allow determining a threshold of dimension of systems equations. Below the threshold the organization of parallel computing adversely effects on the efficiency of computing resources. The paper also studies the effect of a sparse coefficient matrix degree on parallel implementation acceleration.
The given results can be applied not only on computational hydrodynamics and gas dynamics, but also energetics, biology, chemistry, pharmacology, astrophysics, and others.
Keywords: computing aerodynamics, air flow problem, parallel computing, Message Passing Interface MPI.
References
1. Belyaev K.V., Garbaruk A.V., Strelets M.Kh., Shur M.L., Spalart P.R. The experience of turbulence direct numerical simulation on supercomputers. Superkompyuternye dni v Rossii 2016: sb. tr. konf. molod. uchenykh [Proc. Conf. of Young Scientists "Supercomputer days in Russia 2016"]. Moscow, 2016, pp. 357-364. Available at: http://russianscdays.org/files/ pdf16/357.pdf (accessed October 2, 2016).
2. Garbaruk A.V., Strelets M.Kh., Shur M.L. Modelirovanie turbulentnosti v raschetakh slozhnykh techeny [Turbulence Modeling in Complex Flotation Calculations]. St. Petersburg, Technical Univ. Publ., 2012, 88 p.
3. Ortega J. Introduction to Parallel and Vector Solution of Linear Systems. Springer Publ., 1988, 305 p. (Russ. ed.: Moscow, Mir Publ., 1991).
4. Pissanetski S. Sparse Matrix Technology. Academic Press, NY, 1984 (Russ. ed.: Moscow, Mir Publ., 1988, 398 p.).
5. Krukier L.A., Pichugina O.A., Chikina L.G. A review of Krylov subspace methods. XIV Mezhdunar. konf. "Sovremennye problemy matematicheskogo modelirovaniya" [Proc. 14th Int. Conf. "Modern Problems of Mathematical Modeling"]. Rostov-on-Don, SFedU Publ., 2011, 400 p.
6. Balandin M.Yu., Shurina E.P. Metody resheniya SLAU bolshoy razmernosti [Methods of Solving High Dimension SLAE]. Novosibirsk, RGGU Publ., 2000, 70 p.
7. Zhukov V.T., Novikova N.D., Feodoritova O.B. Sdvigovaya strategiya v obobshchennom metode minimalnykh nevyazok [A Shearing Strategy in a Generalized Method of Minimal Residual]. Moscow, 2009, 29 p.
8. Voevodin V.V., Voevodin Vl.V. Parallelnye vychisleniya [Parallel Computing]. St. Petersburg, BHV-Peterburg, 2002.
9. Gergel V.P. Teoriya ipraktikaparallelnykh vychisleny [Theory and Practice ofParallel Computing]. Moscow, BINOM, Laboratoriya znany Publ., 2007, 423 p.
10. Burenkov S.A., Shamayeva O.Yu. An effective function library to solve a class of aerodynamic problems on cluster architectures. XIV Mezhdunar. nauchno-praktich. konf. "Aktualnye voprosy razvitiya innovatsionnoy deyatelnosti v novom tysyacheletii" [Proc. 14th Int. Science and Practice Conf. "Important Issues of Innovation Activities Development in a new Millenium"]. 2015, no. 3 (14), pp. 8-11.
