Software & Systems

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

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

В настоящее время решение задачи аэродинамического обтекания можно получить с использованием методов математического моделирования. Такой подход позволяет с большой степенью точности определить все необходимые характеристики без использования громоздких стендов и проведения дорогостоящих натурных испытаний. Однако моделирование аэродинамических процес- сов относится к классу сложных задач [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)

Получить аналитическое решение системы (1) для большинства практических задач не представляется возможным, поэтому систему Навье–Стокса решают численными методами, для чего вокруг обдуваемого профиля строится расчетная область, изображенная на рисунке 2, которая разбивается неравномерной прямоугольной расчетной сеткой, как показано на рисунке 3.

После дискретизации системы (1) дальнейшее решение задачи сводится к системе линейных алгебраических уравнений (СЛАУ) большой размерности, матрица которой обладает свойствами разреженности и нерегулярности портрета, отсутствием диагонального преобладания и симметрии.

Результаты решения задачи обтекания могут быть представлены графически в виде распределения основных аэродинамических характеристик вокруг обдуваемого профиля. На рисунке 4 изображены распределения давления и скорости в воздушном потоке.

Организация параллельно-последовательных вычислений

На основе анализа методов решения СЛАУ с учетом перечисленных свойств матрицы коэффициентов [3, 4] предпочтение было отдано проекционному методу обобщенных минимальных невязок (General Minimal Residual, GMRES) [5], на основе которого в работе предложена и реализована параллельная модификация GMRES для кластерных архитектур.

Метод обобщенных минимальных невязок относится к классу методов проекционного типа и основывается на проектировании в подпространствах Крылова [5–7].

На рисунке 5 изображена структурная схема параллельного алгоритма решения СЛАУ по методу GMRES, вычисления в которой представлены в виде совокупности взаимодействующих последовательно выполняемых процессов. Для организации взаимодействия используется коммуникационная библиотека передачи сообщений (Message Passing Interface, MPI) [8, 9].

Главный процесс рассылает подчиненным заданные размерность подпространства Крылова m, требуемую точность решения eps, максимальное количество итераций maxNumbIter, а также соответствующие им части матрицы коэффициентов A и вектора правой части b (блоки 2 и 3). В качестве начального приближения к решению выбирается нулевой вектор. Производится предварительная оценка невязки начального приближения r, и вычисляется ее норма betta: i-й процесс производит вычисление фрагмента вектора невязки r(i), после чего фрагменты собираются главным процессом (блок 5).

Далее происходит итерационный процесс вычислений (блоки 6–11), завершающийся при достижении заданной точности eps или при выполнении максимального числа итераций. На каждой итерации процессы совместно формируют базисные вектора w, которые представляются в матричном виде H (блок 8). Полученная матрица H приводится к форме Хессенберга (блок 9) и используется в качестве матрицы коэффициентов треугольной СЛАУ H × y = g (блок 10). Решение полученной треуголь- ной СЛАУ используется при уточнении решения x исходной системы уравнений, полученного на предыдущей итерации. Наконец, по вышеопи- санному принципу вычисляются новая невязка системы r и ее норма betta (блок 11). Отметим, что наиболее вычислительно трудоемкие фрагменты алгоритма выполняются параллельно (бло­ки 5, 8 и 11).

Предложенный алгоритм обладает свойством масштабируемости и предполагает равномерную загрузку всех доступных вычислительных узлов кластера как исходными данными, так и непосредственно вычислительной работой [10].

Оценка трудоемкости параллельного алгоритма

Оценим вычислительную сложность реализованного параллельного алгоритма на множестве вычислительных узлов p при следующих предположениях. Пусть n – размерность решаемой СЛАУ, m – размерность подпространства Крылова, L – количество выполненных итераций, S – длительность аддитивной арифметической операции, P – длительность мультипликативной арифметической операции, d – доля ненулевых элементов в матрице коэффициентов, которая определяется как , где N – количество всех коэффициентов матрицы, а N0 – количество коэффициентов, равных нулю. Тогда трудоемкость T представленного в работе алгоритма может быть оценена с помощью формулы

(2)

В работе выполнено сравнение теоретической оценки трудоемкости алгоритма с экспериментальными данными. На рисунке 6 представлены графики зависимости времени решения СЛАУ с помощью разработанной модификации GMRES от количества используемых вычислительных узлов. Для этой же системы уравнений был проведен вычислительный эксперимент, результаты которого отмечены на графике пунктиром.

Отметим, что время решения СЛАУ в ходе эксперимента немного больше предсказанного оценкой независимо от количества используемых вычислительных узлов. Этот эффект объясняется тем, что оценка (2) не учитывает временные затраты на чтение исходных данных и затраты, необходимые для осуществления коммутационных обменов между исполнителями.

Результаты вычислительного эксперимента

Вычислительный эксперимент проводился на вычислительном кластере из 32 процессоров AMD Opteron 6272 с тактовой частотой 2.1 ГГц, соединенных сетью Gigabit Ethernet, под управлением ОС Linux Kernel 3.12.13. Для оценки эффективности параллельного решения использовалось поня- тие ускорения параллельного решения, определяемого как S(p)=T1/Tp, где T1 – время последовательного решения; Tp – время параллельного решения с использованием p однородных вычислительных узлов. На рисунках 7 и 8 приведены зависимости достигаемых ускорений решений от количества вычислительных узлов. При анализе зависимостей (рис. 7) можно отметить, что для разреженной системы уравнений небольшой размерности (n = 103) ускорение с ростом количества используемых вы- числительных узлов убывает. Это явление закономерно и объясняется тем, что выигрыш от приме- нения распределенных вычислений не компенсирует затраты на их организацию. Для систем уравнений бо́льшей размерности (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.     Жуков В.Т., Новикова Н.Д., Феодоритова О.Б. Сдвиговая стратегия в обобщенном методе минимальных невязок. М.: Изд-во ИПМ РАН, 2009. 28 с.

8.     Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб: БХВ-Петербург, 2002. 608 с.

9.     Гергель В.П. Теория и практика параллельных вычислений. М.: БИНОМ. Лаборатория знаний, 2007. 423 с.

10.   Буренков С.А., Шамаева О.Ю. Библиотека эффективных функций для решения класса аэродинамических задач на кластерных архитектурах // Актуальные вопросы развития инновационной деятельности в новом тысячелетии: сб. тр. XIV Междунар. науч.-практич. конф. 2015. № 3 (14). С. 8–11.