автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Моделирование с помощью МВС двух- и трехмерных течений вязкого газа на основе квазигазодинамических уравнений на нерегулярных сетках

На правах рукописи

аи344 735У Свердлин Александр Александрович

Моделирование с помощью МВС двух- и трехмерных течений вязкого газа на основе квазигазодинамических уравнений на нерегулярных сетках

Специальность 05.13.18 Математическое моделирование, численные методы и комплексы программ

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

з О СЕ И 2505

Москва-2008

003447357

Работа выполнена в Институте математического моделирования Российской академии наук

Научный руководитель:

кандидат физико-математических наук

старший научный сотрудник, зам. директора ИММ РАН

Поляков Сергей Владимирович

Официальные оппоненты:

доктор физико-математических наук, профессор

профессор Российского государственного социального университета

Злотник Александр Анатольевич

кандидат физико-математических наук, доцент старший научный сотрудник ИММ РАН Шильников Евгений Владимирович

Ведущая организация:

Институт прикладной математики им. М.В. Келдыша Российской академии наук

Защита состоится «2-7 'У О/Л 2008 г. в '/7 на заседании

диссертационного совета Д.002.058.01 при Институте математического моделирования РАН по адресу 125047, г. Москва, Миусская пл., 4а.

С диссертацией можно ознакомиться в библиотеке ИММ РАН Автореферат разослан « $I »с^гГ(о 2008 г.

Ученый секретарь диссертационного совета, доктор физико-математических наук, профессор Змитренко Николай Васильевич

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ Актуальность работы

Диссертационная работа посвящена актуальной в настоящее время проблеме создания прикладного программного обеспечения многопроцессорных вычислительных систем (МВС) с распределенной и гибридной памятью (в последнем случае имеется ввиду распределенная система, каждый узел которой представляет собой кластер с общей памятью) для решения современных задач механики сплошной среды, в том числе задач междисциплинарных, сочетающих описания процессов, протекающих в газе, жидкости и твердом теле.

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

Вычислительная гидродинамика (СРБ) зарекомендовала себя лучшим образом для расчета аэродинамического нагрева [1]. Однако особенностями данного семейства методов являются большая вычислительная сложность алгоритмов, большие объемы требуемой памяти, необходимость проведения нового расчета при малейшем изменении параметров течения (например, числа Маха, высоты, угла атаки и т.д.). Поэтому использование этих методов для расчета тепловых потоков на всей траектории полета космического аппарата становится довольно затрудненной и длительной по времени задачей. Перечисленные факторы определяют необходимость использования

высокопроизводительных параллельных систем для моделирования сверхзвуковых течений.

Кроме того, сверхзвуковые течения с большими числами Маха (несколько десятков единиц), характерные для космических полетов, накладывают ограничения на выбор модели течения газа. Предлагаемый в работе квазигазодинамический подход позволяет производить более точные расчеты по сравнению с хорошо известной и широко используемой системой уравнений Навье-Стокса. Также улучшенная вычислительная устойчивость разработанного численного метода позволяет учесть влияние дополнительных факторов, таких как процессы излучения.

Имеющиеся в настоящее время программы расчета вязких течений основаны на использовании системы уравнений Навье-Стокса. Предлагаемая в работе математическая модель использует квазигазодинамические (КГД) уравнения [2,3,4], которые отличаются от уравнений Навье-Стокса дополнительными диссипативными слагаемыми с малым параметром в качестве коэффициента. Использование дополнительной диссипации позволяет существенно оптимизировать вычислительные алгоритмы. Это свойство делает особенно выигрышными численные алгоритмы решения задач радиационной газовой динамики, основанные на КГД уравнениях [5]. Это связано с тем, что радиационные добавки в уравнении энергии вносят мощный дополнительный фактор неустойчивости. В системе уравнений Навье-Стокса не имеется возможностей для эффективной регуляризации численных решений в этом случае. Система КГД уравнений, напротив, обладает естественным регуляризатором.

Работа выгодно отличается от имеющихся в свободном доступе программных средств по следующим критериям:

• возможность задания реальной геометрии объекта с произвольной точностью;

• использование неструктурированных расчетных сеток (двухмерные гибридные и треугольные, трехмерные блочные и тетраэдральные);

• использование системы квазигазодинамических уравнений в качестве математической модели;

• расчет течений в широком диапазоне скоростей потока - от дозвукового, до сверхзвукового (50Ма и более);

• возможность учета дополнительных факторов (например, учет процессов излучения в газе, химические реакции, многокомпонентная среда - газ + жидкость + твердое тело);

• эффективное использование МВС с распределенной и гибридной архитектурой различной производительности - от бюджетных кластеров (32-64 процессора, объединенных простой сетью Ethernet) до высокопроизводительных параллельных систем (2000-4000 и более процессоров, объединенных сетями Myrinet, Infiniband).

Цели и задачи диссертационной работы

• Разработать численный метод на основе системы КГД уравнений для моделирования сверхзвуковых двух- и трехмерных течений вязкого газа вокруг тел сложной формы на нерегулярных сетках большого объема.

• Разработать высокоэффективный параллельный алгоритм решения задачи и реализовать его в виде программного комплекса для высокопроизводительных систем с распределенной и гибридной памятью.

• Провести расчеты двух- и трехмерных модельных задач, сравнить полученные результаты с известными данными.

Научная новизна и практическая ценность работы

В диссертации предложены оригинальные численные методы решения задач газовой динамики на основе аппроксимации КГД уравнений на двухмерных гибридных (на основе прямоугольных и треугольных ячеек) и трехмерных тетраэдральных сетках. Предложен и реализован оригинальный алгоритм генерации адаптивных гибридных локально-сгущающихся расчетных сеток в двухмерной области. Предложены и реализованы

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

Проведено численное моделирование ряда модельных задач. Полученные данные сопоставлены с результатами, полученными другими авторами. Качественное совпадение с теоретическими расчетами [2] и экспериментальными данными [6,7] подтверждает достоверность полученных результатов. Реализация результатов

Результаты диссертации получены в ходе выполнения работ по проектам РФФИ, в которых автор принимал участие в качестве исполнителя (гранты № 05-07-90230, 06-01-00233,08-07-00458). Личный вклад автора

Автором работы получены аппроксимации системы КГД уравнений на неструктурированных двухмерных гибридных и трехмерных тетраэдральных сетках. Разработаны эффективные параллельные алгоритмы на основе явных

6

численных схем, использующих данные аппроксимации. Выполнены параллельные реализации алгоритмов с использованием библиотеки MPI для архитектур с распределенной памятью и технологии ОрепМР - для систем с общей памятью. Созданы генераторы локально-сгущающихся адаптивных неструктурированных сеток. Для двухмерных областей предложены собственные оригинальные алгоритмы. Для генерации трехмерных тетраэдральных сеток используется библиотека TetGen, однако, локальное сгущение и адаптация сетки к особенностям решения осуществляется с помощью алгоритмов, разработанных автором диссертации. Для задания геометрии исследуемых объектов предложена технология представления объекта в виде сплайнов и кусочно-линейных элементов, реализованная в виде пр01раммы - редактора геометрии. Апробация работы

Основные результаты диссертационной работы докладывались и обсуждались на многих российских и международных научно-технических конференциях и семинарах, в частности:

1. Международная конференция Parallel CFD 2005, 24-27 мая 2005 г., University of Maryland, США. Устный доклад «Numerical Simulation of 2D Radiative Heat Transfer for Reentry Vehicles» (соавторы C.B. Поляков, T.A. Кудряшова, Э.М. Кононов)

2. Международная конференция «Тихонов и современная математика», 19-25 июня 2006 г., Москва, МГУ им. М.В. Ломоносова. Устный доклад «Численное моделирование двумерных задач переноса радиации» (соавторы C.B. Поляков, Т.А. Кудряшова, Э.М. Кононов).

3. Международная конференция «Тихонов и современная математика», 19-25 июня 2006 г., Москва, МГУ им. М.В. Ломоносова. Устный доклад «Расчет параметров течения газа вокруг спускаемого аппарата».

4. Международная конференция PARA'06, 18-21 июня 2006 г., Umea University, Швеция. Устный доклад «Numerical Simulation of 2D

Radiation Heat Transfer Problem» (соавторы C.B. Поляков, Т.А. Кудряшова, Э.М. Кононов).

5. Международная конференция WEHSFF 2007, 19-22 ноября 2007 г., Москва. Устный доклад «Parallel Computation of Radiation Transport around Reentry Vehicle» (соавторы C.B. Поляков, Т.А. Кудряшова, Э.М. Кононов).

6. Пятый международный научный семинар «Математические модели и моделирование в лазеро-плазменных процессах», 29 января - 2 февраля 2008 г., Кусково, МО. Устный доклад «Численное моделирование проблемы переноса излучения вокруг возвращаемого аппарата» (соавторы C.B. Поляков, Т.А. Кудряшова, Э.М. Кононов).

7. Международная конференция Parallel CFD 2008, 19-22 мая 2008 г., Lyon, Франция. Устный доклад «3D Numerical Simulation of Gas Flow Around Reentry Vehicles» (соавторы C.B. Поляков, T.A. Кудряшова, Э.М. Кононов).

Публикации

Основные результаты диссертации опубликованы в трех рецензируемых работах (см. список в конце автореферата) и тезисах выше указанных конференций.

Структура и объем работы

Диссертация состоит из введения, четырех глав, заключения и списка литературы. Объем составляет 110 машинописных страниц, текст содержит 38 рисунков и 16 таблиц.

СОДЕРЖАНИЕ РАБОТЫ Во введении обоснована актуальность рассматриваемых в работе проблем, сформулированы основные цели диссертации и дано ее краткое содержание по главам. Показана важность вычислительного эксперимента в аэродинамике. Приведен краткий обзор существующих методов и подходов к моделированию вязких сжимаемых течений. Дана краткая классификация

современных многопроцессорных вычислительных систем. Рассмотрены различные подходы к параллелизации вычислительных алгоритмов.

Первая глава посвящена описанию решаемой задачи и соответствующей ей математической модели для случаев двухмерного и трехмерного течений.

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

Рисунок 1. Модель спускаемого аппарата Apollo.

Параметры течения газа определяются на основе известных характеристик атмосферы Земли на высоте полета, согласно траектории спускаемого аппарата. Параметры течений, для которых проводились исследования, приведены в четвертой главе. Проведено численное моделирование ряда задач, для которых доступны экспериментальные данные [6,7]. В работе рассматривались задачи двухмерного осесимметричного и трехмерного обтекания летательного аппарата. В случае двухмерной задачи, использовалась прямоугольная расчетная область. Для моделирования трехмерных задач задавалась замкнутая область произвольной формы, границы которой представимы в виде совокупности треугольников в пространстве XYZ.

Для описания течения газа использовалась система (1)-(3) квазигазодинамических (КГД) уравнений [2,3,4], дополненная граничными

9

условиями, уравнениями состояния газа, выражениями для коэффициентов вязкости (4) и теплопроводности (5). КГД уравнения обобщают хорошо известную систему уравнений Навье-Стокса (НС) и отличаются от нее дополнительными диссипативными слагаемыми, с мальм параметром т. Эти слагаемые могут бьггь охарактеризованы как искусственная вязкость и обеспечивают дополнительную вычислительную устойчивость алгоритма.

1+^=0 (1) е(ри).

з/

- + ® и)+ Ур = сНуП

и р 2 р

(2)

+ = (3)

= —[с!Ь'(ри ® и)+ Vр] Р

П = ПЮ + та ® [¿>(и • У)и + Ур]+т/[(и • У)р + десЦуи] (4)

4 = %да-три

(5)

КГД уравнения связывают между собой газодинамические величины плотности, импульса и энергии. Дополнительные величины температуры и давления вычисляются в соответствии с уравнениями состояния газа.

На части границы расчетной области задаются условия втекающего (невозмущенного) потока. В случае моделирования задач, имеющих явную симметрию, на одной или нескольких частях границы задаются соответствующие граничные условия. На оставшейся части границы задаются «мягкие» условия, позволяющие потоку газа свободно вытекать из расчетной области. Граничные условия на поверхности объекта не позволяют потоку газа проникать внутрь. Течение считается квазистационарным, расчет проводится на установление решения. В качестве начальных условий используются параметры невозмущенного потока [2].

Во второй главе рассматривается численный алгоритм расчета сверхзвуковых течений.

Конечно-разностные аппроксимации КГД системы строятся в потоковой форме, используя вектор плотности потока массы .¡щ, тензор вязких напряжений П и вектор теплового потока q, соответствующие законам сохранения (1)-(3). Инвариантная форма записи КГД системы позволяет конструировать численные схемы для любой ортогональной системы координат для структурированных и неструктурированных расчетных сеток.

Для численного решения задачи используются сетки по времени и пространству. Газодинамические параметры (р, и, Е, р, Т) определяются в узлах пространственной сетки. Значения газодинамических переменных в центрах ячеек определяются как арифметические средние значений в прилегающих узлах.

Для пространственной дискретизации используются неструктурированные адаптивные локально-сгущающиеся к особенностям решения сетки (треугольные и гибридные в двухмерном случае и тетраэдральные в трехмерном). Локальное сгущение позволяет оптимально сбалансировать точность расчета и вычислительную сложность и является ключевым фактором для трехмерного расчета, характеризующегося большим числом узлов сетки даже для простейших задач. При расчете модельных задач использовались сетки, содержащие от 100 тысяч до 2 миллионов узлов. Процедура подготовки сетки для параллельных расчетов подробно описана в третьей главе.

Явная по времени численная схема, основанная на КГД уравнениях, строится методом контрольных объемов (МКО). Значения газодинамических переменных определяются в узлах сетки, вокруг которых строятся барицентрические контрольные объемы (Рисунки 2,3,4). Контрольные объемы полностью покрывают всю расчетную область без промежутков и пересечений. Пространственные производные аппроксимируются со вторым порядком точности, производные по времени - с первым порядком точности.

/V

Рисунок 3. Тетраэдральная ячейка и грани барицентрического контрольного объема

Рисунок 4. Регион тетраэдральной сетки и пример контрольного объема Ограничение на временной шаг для такого рода схем определяется условием Куранта, которое имеет вид:

Д? = А™ (6)

^шах

где - минимальное ребро треугольника в сетке, (3 - числовой

коэффициент. Критерий завершения расчета выбран на основе величины приращения плотности на временном шаге:

\Pnexl Р рпч

!-'-<В {/)

Д/

В третьей главе подробно освещаются алгоритмы параллелизации, описывается состав и функции программного комплекса.

Параллелизация КГД вычислительного кода производится с использованием метода геометрического параллелизма. Вся расчетная область разбивается на подобласти с примерно равным количеством узлов сетки. Количество таких подобластей соответствует количеству вычислительных ядер конкретной параллельной системы, которые планируется использовать в данном запуске. Разбиение на подобласти проводится таким образом, чтобы минимизировать количество общих (граничных) узлов между подобластями и, как следствие, объем передаваемой между вычислителями информации. Далее подобласти распределяются между различными узлами кластера в случае систем с распределенной памятью, а так же между вычислительными ядрами системы с общей памятью, причем для гибридных систем (когда один вычислительный узел распределенного кластера является параллельной системой с общей памятью) применяются оба типа распараллеливания. Для распараллеливания в рамках модели распределенной памяти применяется библиотека MPI, для модели общей памяти используется библиотека ОрепМР. Все расчетные модули программного комплекса написаны на С++ и являются кросс-платформенными.

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

Считывание начальных данных для расчета осуществляется на главном (нулевом) узле вычислительной сети, после чего данные передаются на остальные вычислительные модули в соответствии с заранее определенным разбиением расчетной области. Периодически, с заранее определенным интервалом, либо при выполнении критерия завершения расчета

выполняется сбор информации со всех узлов и сохранение ее в файлах результатов на главном узле вычислительной системы.

Для измерения эффективности параллельных алгоритмов использовалась следующая методика. Производился замер среднего времени, затрачиваемого на выполнение одного временного шага по КГД уравнениям. Усреднение проводилось на интервале в 4096 шагов. Такое среднее «шаговое» время на исследуемой конфигурации (количестве процессоров) сопоставлялось со средним «шаговым» временем при использовании только одного узла вычислительной системы (для систем МВС-ЮОК и СКИФ один узел - это 8 процессоров, для системы IMM2 - два процессора). Для расчета эффективности распараллеливания применялась следующая формула:

Efficiency = -L-~ '100% , (8) t„-N

где tj - среднее время итерации при использовании одного узла ВС, Р -количество процессоров на одном узле; tN - среднее время итерации на N процессорах, N - количество процессоров, на которых запускается исследуемая конфшурация.

Характеристики использованных вычислительных систем: Кластер IMM2 (ИММ РАН): 26 вычислительных узлов, включающих по 2 процессора Intel Xeon ЗГТц; пиковая производительность 300 GFlop/s; пиковая пропускная способность 100 Мбит/с.

Суперкомпьютер МВС-ЮОК (МСЦ РАН): 470 вычислительных узлов, включающих по 2 процессора Intel Xeon ЗГТц (8 ядер) и 4Гб ОЗУ; пиковая производительность 45 TFlop/s; пиковая пропускная способность 1540 Мбайт/с.

Суперкомпьютер СКИФ (МГУ): 630 вычислительных узлов, включающих по 2 процессора Intel Xeon ЗГТц (8 ядер) и 8Гб ОЗУ; пиковая производительность 60 TFlop/s; пиковая пропускная способность 1540 Мбайт/с.

Эффективность параллельного кода исследовалась на системах различной производительности, от бюджетных кластеров, состоящих из нескольких десятков процессоров, объединенных сетью Ethernet, до высокопроизводительных вычислительных систем на базе нескольких тысяч процессоров, объединенных сетями Myrinet, Infiniband. Программный комплекс продемонстрировал высокую эффективность параллельных алгоритмов, особенно для задач большого объема. Ниже приводятся графики эффективности параллелизации КГД алгоритма на различных параллельных вычислительных системах (Рисунки 5,6,7):

Число процессоров

Рисунок 5. Анализ эффективности параллелизации КГД алгоритма для 20 задач на вычислительном комплексе МВС-ЮОК (МСЦ РАН)

Число процессоров

Рисунок 6. Анализ эффективности параллелизации КГД алгоритма для ЗБ задач средствами МИ на вычислительном комплексе СКИФ (МГУ)

Эффективность параллелизации КГД, сетка 1.1 M узлов

Кол-во процессоров

Рисунок 7. Сравнительная эффективность параллелизации КГД алгоритма для 3D задач средствами MPI и MPI+OpenMP на вычислительном комплексе МВС-ЮОК (МСЦ РАН)

16

В конце третьей главы подробно описывается состав программного комплекса, излагаются детали реализации алгоритмов, выполнен обзор программных средств обработки и визуализации результатов расчета. Разработанный программный комплекс включает:

• Редактор геометрии (для задания расчетной области и исследуемого объекта)

• Программу генерации адаптивных локально сгущающихся неструктурированных сеток

• Расчетный модуль для параллельных ВС, включающий блок расчета газодинамических уравнений и блок расчета радиационной газовой динамики

Основными форматами файлов результатов являются файл .node (библиотеки TetGen) и .pit файлы для визуализации в TecPlot.

В четвертой главе приводятся результаты моделирования тестовых задач. Проведено сопоставление с известными результатами точных решений и экспериментов. Приведем параметры одной из модельных задач:

- Температура воздуха в невозмущенном потоке (Т0): 266 К

- Температура внутри фронта ударной волны (Т(): 2500-5500 К

- Давление воздуха в невозмущенном потоке (Р0): 43 Па

- Плотность воздуха в невозмущенном потоке (р0): 5,63-Ю"4 кг/м3

- Характерный размер задачи (Lo): 4 м -Число Маха: 12

Пример расчета показан на рисунке 8, изображены изолинии числа Маха вокруг объекта.

z

Рисунок 8. Изолинии числа Маха вокруг модели спускаемого аппарата Apollo при скорости набегающего потока 12 Ма В заключении приведены основные результаты диссертации, выносимые на защиту; сформулированы основные выводы и обозначены перспективы для дальнейшей работы.

ЛИТЕРАТУРА

1. Gnoffo Р.А., K.J. Weilmuenster, H.H.Hamilton II, D.R. Olynick, E. Venkatapathy. Computational Aerothermodynamic Design Issues for Hypersonic Vehicles, Journal of Spacecraft and Rockets, vol. 36, no. 1, January-February 1999, pp. 21^43.

2. Т.Г. Елизарова. Квазигазодинамические уравнения и методы расчета вязких течений. - М: Научный мир, 2007. - 352 с.

3. Шеретов Ю.В. Анализ устойчивости модифицированной кинетически-согласованной схемы в акустическом приближении // Применение функционального анализа в теории приближений. - Тверь: Тверской государственный университет, 2004. - 147-160.

4. Елизарова Т.Г., Соколова М.Е. Численный алгоритм расчета сверхзвуковых течений, основанный на квазигазодинамических уравнениях. Вестник Московского университета, серия 3. Физика. Астрономия, 2004, N 1, с. 10-15.

5. Т.А. Кудряшова, C.B. Поляков, Э.М. Кононов. Расчёт поля радиационного излучения газа вокруг спускаемого аппарата. // Математическое моделирование, 2008, 20(10), с. 1-12 (принята к печати).

6. William L. Kleb, К. James Weilmuenster. A High Angle of Attack Inviscid Shuttle Orbiter Computation. NASA TM 107606, Apr. 1992.

7. William L. Ко, Leslie Gong, Robert D. Quinn. Hypothetical Reentry Thermostructural Performance of Space Shuttle Orbiter With Missing or Eroded Thermal Protection Tiles. NASA TM-2004-212850. July 2004.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

1. Предложен новый параллельный численный алгоритм решения газодинамической задачи на основе КГД уравнений на нерегулярных сетках большой размерности.

2. Предложена методика задания областей сложной формы с помощью сплайнов.

3. Разработаны алгоритмы генерации сеток различного типа: треугольных, тетраэдральных и гибридных (сочетающих ячейки треугольной и прямоугольной формы).

4. Разработанные алгоритмы реализованы в виде комплекса параллельных программ для МВС различной архитектуры.

5. Получены результаты расчетов модельных двух- и трехмерных задач на сетках различного типа.

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

1. B.N. Chetverushkin, S.V. Polyakov, Т.А. Kudryashova, Е.М. Kononov, АА. Sverdlin, Numerical Simulation of 2D Radiative Heat Transfer for Reentry Vehicles, Proc. Parallel CFD 2005 Conference, Elsevier, 2006, pp. 293-299.

2. Sergey Polyakov, Tatiana Kudryashova, Alexander Sverdlin, Eldar Kononov. Parallel Computation of Radiation Transport around Reentry Vehicle. / CD-proceedings of "West-East High Speed Flow Field Conference (WEHSFF 2007)" (November 19-22, 2007, Moscow, Russia), 2007, pp. 1-8.

3. Т.А. Кудряшова, C.B. Поляков, A.A. Свердлин. Расчёт параметров течения газа вокруг спускаемого аппарата. // Математическое моделирование, 2008,20(8), с. 119-128.

Напечатано с готового оригинал-макета

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. Подписано к печати 04.09.2008 г. Формат 60x90 1/16. Уел печ.л. 1,25 Тираж 100 экз. Заказ 482. Тел. 939-3890. Тел./факс 939-3891. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. МБ. Ломоносова, 2-й учебный корпус, 627 к.

Введение

1 КГД модель сверхзвуковых течений газа и постановки задач

1.1 Постановка модельной задачи.

1.1.1 Исследуемый объект.

1.1.2 Параметры течения

1.2 Система КГД уравнений.

1.2.1 Обезразмеривание КГД системы.

1.2.2 Введение искусственной диссипации.

1.3 Моделирование двухмерных течений.

1.3.1 Геометрия расчетной области.

1.3.2 Система КГД уравнений в RZ геометрии.

1.3.3 Начальные и граничные условия для двухмерных задач.

1.4 Моделирование трехмерных течений.

1.4.1 Геометрия расчетной области.

1.4.2 Система КГД уравнений в XYZ геометрии.

1.4.3 Начальные и граничные условия для трехмерных расчетов.

1.5 Выводы.

2 Численный алгоритм расчета сверхзвуковых течений

2.1 Численный алгоритм расчета двухмерных течений.

2.1.1 Расчетная область и сетка.

2.1.2 Разностная аппроксимация уравнений.

2.1.3 Аппроксимация частных производных.

2.1.4 Аппроксимация начальных и граничных условий

2.1.5 Алгоритм расчета

2.1.6 Устойчивость и точность численной схемы.41"

2.2 Численный алгоритм расчета трехмерных течений.

2.2.1 Расчетная область и сетка.

2.2.2 Разностная аппроксимация уравнений.

2.2.3 Аппроксимация частных производных.

2.2.4 Аппроксимация начальных и граничных условий

2.2.5 Критерий установления течения.

2.2.6 Устойчивость и точность численной схемы.

2.3 Выводы.

3 Алгоритмы распараллеливания и программный комплекс

3.1 Параллельный алгоритм для систем с распределенной памятью

3.1.1 Подготовка сетки для параллельных расчетов и балансировка загрузки.

3.1.2 Построение схемы пересылки данных.

3.1.3 Обмен данными (библиотека MPI).

3.1.4 Параллельный ввод-вывод данных.

3.2 Описание оборудования.

3.3 Анализ эффективности параллельного алгоритма.

3.3.1 Метод измерения эффективности.

3.3.2 Эффективность параллельного алгоритма

3.3.3 Анализ результатов

3.4 Оптимизация параллельного алгоритма для систем с гибридной архитектурой

3.4.1 Технология ОрепМР.

3.4.2 Гибридный алгоритм параллелизации.

3.4.3 Эффективность гибридного алгоритма

3.5 Комплекс программ.

3.5.1 Программа редактирования двухмерной геометрии

3.5.2 Программа генерации двухмерных расчетных сеток

3.5.3 Программа моделирования двухмерных задач обтекания

3.5.4 Программа генерации тетраэдральных сеток

3.5.5 Программа моделирования трехмерных задач обтекания

3.5.6 Визуализация трехмерной расчетной области и сетки: TetView.

3.5.7 Визуализация результатов расчета: Tecplot.

4 Результаты моделирования

4.1 Моделирование двухмерного течения

4.2 Моделирование трехмерных течений.

4.2.1 Модель спускаемого аппарата Apollo.

4.2.2 Орбитальный корабль Space Shuttle.

Диссертационная работа посвящена актуальной в настоящее время проблеме создания прикладного программного обеспечения многопроцессорных вычислительных систем (МВС) с распределенной и гибридной памятью (в последнем случае имеется ввиду распределенная система, каждый узел которой представляет собой систему с общей памятью) для решения современных задач механики сплошной среды, в том числе задач междисциплинарных, сочетающих описания процессов, протекающих в газе, жидкости и твердом теле. В качестве основной темы в работе выбрана проблема моделирования сверхзвуковых течений вязкого газа на сетках различной структуры. В качестве основного приложения исследуется задача о расчете параметров течения газа вблизи обшивки спускаемого в атмосфере Земли космического аппарата (КА). Начнем введение в проблему с указанного выше приложения.

Численное моделирование аэротермодинамических процессов играет огромную роль в конструировании любого космического корабля и, в частности, спускаемого аппарата. При его проектировании крайне важно знать распределение температуры на поверхности КА, чтобы определить тепловые нагрузки и убедиться, что они не выходят за рамки конструктивных ограничений и не приведут к его разрушению. Тепловая нагрузка на структуру летательного аппарата должна учитываться наравне с аэродинамическими нагрузками. Начиная с первых космических запусков, разработка теплозащиты была "Проблемой №1" [1]. Первоначально масса теплозащиты доходила до 30% массы космического корабля. Последние эксперименты показывают, что тепловая защита космических аппаратов, разработанная десятилетия назад, оказывается избыточной, и величины разрушительных факторов, таких как температура, не достигают тех значений, на которые тепловая защита была рассчитана [2]. Программы выведения коммерческих космических аппаратов предъявляют более жесткие требования к расчету газодинамических величин в полете. Моделирование [3] позволяет экономить огромные средства, а также предсказывать параметры газодинамических процессов, которые невозможно получить в условиях экспериментов на Земле.

Методы вычислительной гидродинамики (CFD ) зарекомендовали себя лучшим образом для расчета аэродинамического нагрева [4]. Однако особенностями данного семейства методов являются большая вычислительная сложность алгоритмов, большие объемы требуемой памяти, необходимость проведения нового расчета при малейшем изменении конструкции КА и параметров течения (например, числа Маха, высоты, угла атаки и т.д.). Поэтому использование этих методов для расчета тепловых потоков на всей траектории полета КА затрудняется большим временем вычислений. Перечисленные факторы определяют необходимость использования высокопроизводительных параллельных систем для моделирования сверхзвуковых течений.

Кроме того, сверхзвуковые течения с большими числами Маха (несколько десятков единиц), характерные для космических полетов, накладывают ограничения на выбор модели течения газа. Имеющиеся в настоящее время программы расчета вязких течений основаны на использовании системы уравнений Навье-Стокса. Предлагаемая в работе математическая модель использует квазигазодинамические (КГД) уравнения [5, 6, 7], которые отличаются от уравнений Навье-Стокса дополнительными диссипативными слагаемыми с малым параметром в качестве коэффициента. Это свойство делает особенно предпочтительными численные алгоритмы решения задач радиационной газовой динамики (РГД), основанные на КГД уравнениях [8]. Последнее связано с тем, что радиационные слагаемые в уравнении энергии вносят мощный дополнительный фактор неустойчивости. В системе уравнений Навье-Стокса не имеется возможностей для эффективной регуляризации численных решений в таких случаях. Система КГД уравнений, напротив, обладает естественным регуляризатором. В случае стационарных течений КГД слагаемые практически не влияют на точность решения, но, как показано в [5], обеспечивают более точное решение для течений с ударными волнами и резкими градиентами газодинамических величин.

Моделированию течений газа на основе КГД уравнений посвящены работы А.Е. Дуйсекулова, И.А. Граур [9, 10], И.А. Широкова, М.Е. Соколовой [11, 12, 13], связанные с расчетами на регулярных сетках. Расчеты с использованием КГД системы на треугольных сетках проводились в работах В.В. Серегина [14]. Работы, связанные с расчетами по родственным кинетически согласованным разностным схемам (КСРС) [15] на регулярных сетках, проводили И.А. Граур [16, 17, 18], А.Н. Антонов [19], Е.В. Шильников [20, 21], И.В. Абалакин и А.В. Жохова [22], М.А. Шумков [23, 24]. Расчеты на треугольных сетках с использованием КСРС присутсвуют в работах А.В. Жоховой [25] и И.В. Абалакина [26, 27], на тетраэдральных сетках - в работах С.А. Сукова [28] и И.В. Абалакина [29]. Метод конечных элементов применительно к газодинамическим расчетам на тетраэдральных сетках используется в работах С.А. Сукова, И.В. Абалакина, А.В. Горобца [30]. Из анализа этих работ становится ясно, что систематические расчеты на основе КГД уравнений на нерегулярных сетках произвольной структуры не проводились. Нерешенные проблемы существовали при использовании треугольных сеток в случае (r,z) геометрии и при использовании тетраэдральных сеток в трехмерном случае. Определенные сложности были связаны с распараллеливанием таких численных подходов.

В связи со сказанным целями данной диссертационной работы были адаптация КГД подхода к вычислениям на нерегулярных сетках и создание программного комплекса для моделирования двух- и трехмерных течений газа, предназначенного для высокопроизводительных вычислительных систем. Для достижения этих целей были поставлены и успешно решены следующие задачи:

• Разработка численного метода на основе системы КГД уравнений для моделирования сверхзвуковых двухмерных и трехмерных течений вязкого газа вокруг тел сложной формы на нерегулярных сетках большого объема.

• Разработка высокоэффективного параллельного алгоритма решения выбранного класса задач и реализация его в виде программного комплекса для высокопроизводительных систем с распределенной и гибридной памятью.

• Проведение расчетов двухмерных и трехмерных модельных задач, сравнение полученных результатов с известными данными.

Разработанный в рамках диссертационной работы программный комплекс отличается от имеющихся в свободном доступе программных средств рядом преимуществ, а именно:

• возможность задания реальной геометрии объекта с произвольной точностью;

• использование неструктурированных расчетных сеток (двухмерные гибридные и треугольные, трехмерные блочные и тетраэдральные);

• использование системы квазигазодинамических уравнений в качестве математической модели;

• расчет течений в широком диапазоне скоростей потока - от дозвукового до сверхзвукового (50 Ма и более);

• возможность учета дополнительных факторов (например, учет процессов излучения в газе, химических реакций, многокомпонентная среда - газ, жидкость, твердое тело);

• эффективное использование МВС с распределенной и гибридной архитектурой различной производительности - от бюджетных кластеров (32-64 процессора, объединенных простой сетью Ethernet) до высокопроизводительных параллельных систем (2000-4000 и более процессоров, объединенных сетями Myrinet, Infiniband).

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

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

Для проверки корректности численного алгоритма проведено численное моделирование ряда модельных задач. Полученные данные сопоставлены с результатами, полученными другими авторами. Качественное совпадение с теоретическими расчетами [5] и экспериментальными данными [31] подтверждает достоверность полученных результатов.

В рамках диссертационной работы получены аппроксимации системы КГД уравнений на неструктурированных двухмерных гибридных и трехмерных тетраэдральных сетках. Разработаны эффективные параллельные алгоритмы на основе явных численных схем, использующих данные аппроксимации. Выполнены параллельные реализации алгоритмов с использованием библиотеки MPI [32, 33] для архитектур с распределенной памятью и технологии ОрепМР [34] - для систем с общей памятью. Созданы генераторы локально-сгущающихся адаптивных неструктурированных сеток. Для двухмерных областей предложены собственные оригинальные алгоритмы. Для генерации трехмерных тетраэдральных сеток используется библиотека TetGen [35], однако, локальное сгущение и адаптация сетки к особенностям решения осуществляется с помощью алгоритмов, разработанных специально для решения поставленных в диссертации задач. Для задания геометрии исследуемых объектов предложена технология представления объекта в виде сплайнов и кусочно-линейных элементов, реализованная в виде программы - редактора геометрии.

Основные положения, выносимые на защиту:

1. Предложен новый параллельный численный алгоритм решения газодинамической задачи на основе КГД уравнений на нерегулярных сетках большой размерности.

2. Предложена методика задания областей сложной формы с помощью сплайнов.

3. Разработаны алгоритмы генерации сеток различного типа: треугольных, тетраэдральных и гибридных (сочетающих ячейки треугольной и прямоугольной формы).

4. Разработанные алгоритмы реализованы в виде комплекса параллельных программ для МВС различной архитектуры.

5. Получены результаты расчетов модельных двухмерных и трехмерных задач на сетках различного типа.

Основные результаты, представленные в диссертационной работе, докладывались и обсуждались на многих российских и международных научно-технических конференциях и семинарах, в частности:

1. Международная конференция Parallel CFD 2005, 24-27 мая 2005 г., University of Maryland, США. Устный доклад "Numerical Simulation of 2D Radiative Heat Transfer for Reentry Vehicles" (соавторы С.В. Поляков, Т.А. Кудряшова, Э.М. Кононов).

2. Международная конференция "Тихонов и современная математика", 19-25 июня 2006 г., Москва, МГУ им. М.В. Ломоносова. Устный доклад "Численное моделирование двумерных задач переноса радиации" (соавторы С.В. Поляков, Т.А. Кудряшова, Э.М. Кононов).

3. Международная конференция "Тихонов и современная математика", 19-25 июня 2006 г., Москва, МГУ им. М.В. Ломоносова. Устный доклад "Расчет параметров течения газа вокруг спускаемого аппарата".

4. Международная конференция PARA'06, 18-21 июня 2006 г., Umea University, Швеция. Устный доклад "Numerical Simulation of 2D Radiation Heat Transfer Problem" (соавторы С.В. Поляков, Т.А. Куд-ряшова, Э.М. Кононов).

5. Международная конференция WEHSFF 2007, 19-22 ноября 2007 г., Москва. Устный доклад "Parallel Computation of Radiation Transport around Reentry Vehicle" (соавторы С.В. Поляков, Т.А. Кудря-шова, Э.М. Кононов).

6. Пятый международный научный семинар "Математические модели и моделирование в лазеро-плазменных процессах", 29 января - 2 февраля 2008 г., Кусково, МО. Устный доклад "Численное моделирование проблемы переноса излучения вокруг возвращаемого аппарата" (соавторы С.В. Поляков, Т.А. Кудряшова, Э.М. Кононов).

7. Международная конференция Parallel CFD 2008, 19-22 мая 2008 г., Лион, Франция. Устный доклад "3D Numerical Simulation of Gas Flow Around Reentry Vehicles" (соавторы С.В. Поляков, Т.А. Кудряшова, Э.М. Кононов).

8. Международная конференция ECCOMAS 2008, 30 июня - 5 июля 2008 г., Венеция, Италия. Устный доклад " Numerical Simulation of Radiation Transport Problem around Reentry Vehicle" (соавторы С.В. Поляков, Т.А. Кудряшова, Э.М. Кононов).

Результаты, представленные в диссертационной работе, опубликованы в тезисах вышеуказанных конференций и в трех рецензируемых работах:

1. B.N. Chetverushkin, S.V. Polyakov, Т.А. Kudryashova, Е.М. Kononov, А.А. Sverdlin, Numerical Simulation of 2D Radiative Heat Transfer for Reentry Vehicles, Proc. Parallel CFD 2005 Conference, Elsevier, 2006, pp. 293-299.

2. Sergey Polyakov, Tatiana Kudryashova, Alexander Sverdlin, Eldar Kononov. Parallel Computation of Radiation Transport around Reentry Vehicle. / CD-proceedings of "West-East High Speed Flow Field

Conference (WEHSFF 2007)" (November 19-22, 2007, Moscow, Russia), 2007, pp. 1-8.

3. T.A. Кудряшова, С.В. Поляков, А.А. Свердлин. Расчёт параметров течения газа вокруг спускаемого аппарата. // Математическое моделирование, 2008, 20(8), с. 119-128.

Диссертация состоит из введения, четырех глав, заключения и списка литературы.

2.3 Выводы

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

Глава 3

Алгоритмы распараллеливания и программный комплекс

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

Адаптация алгоритма моделирования течения газа на основе КГД уравнений для параллельных систем с распределенной памятью осуществляется по принципу геометрического параллелизма. Обмен данными между отдельными вычислительными узлами осуществляется средствами библиотеки MPI (Message Passing Interface), которая стала стандартом де-факто для многопроцессорных параллельных систем любого уровня, построенных на различных платформах. Подход к параллелиза-ции алгоритмов расчета двухмерных и трехмерных течений одинаков и различается только низкоуровневыми структурами данных.

3.1 Параллельный алгоритм для систем с распределенной памятью

Метод декомпозиции расчетной области [51] предназначен для решения краевой задачи путем разбиения ее на меньшие краевые задачи на подобластях с последующим разрешением границы. Внутри каждой подобласти задача решается независимо, что делает методы декомпозиции расчетной области подходящими для параллельного программирования.

Параллелизация алгоритмов решения краевых задач, основанная на методе декомпозиции расчетной области, называется методом геометрического параллелизма.

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

Упрощенная блок-схема исходного последовательного алгоритма расчета приведена на рис. 3.1(a). Дополнительно в программу должны быть добавлены следующие функции:

• передача считанных из файла данных о сетке и состоянии газодинамических переменных на все вычислительные узлы;

• разбиение расчетной сетки на подобласти;

• построение графа межпроцессорных обменов;

• передача значений газодинамических величин в граничных узлах между вычислителями, обрабатывающими смежные подобласти;

• сбор результирующих данных со всех вычислителей для сохранения результатов работы в файлах.

Упрощенная блок-схема программы, осуществляющей параллельные расчеты, приведена на рис. 3.1(6).

3.1.1 Подготовка сетки для параллельных расчетов и балансировка загрузки

Пространственная расчетная сетка для моделирования газодинамических процессов представляет собой набор узлов (в которых задаются

НАЧАЛО

КОРНЕВОЕ ЯДРО? MPICommrank()

СЧИТЫВАНИЕ ФАЙЛОВ: УЗЛОВ (.node) ТЕТРАЭДРОВ (.elel/

ПЕРЕДАЧА I ПРИЕМ УЗЛОВ I ТЕТРАЭДРОВ . MPLBcastQ ^

РАЗБИЕНИЕ НА ПОДОБЛАСТИ, ПОСТРОЕНИЕ ГРАФА ОБМЕНОВ

ИТЕРАЦИЯ КГД АЛГОРИТМА

ПЕРЕДАЧА /ПРИЕМ ГРАНИЧНЫХ УЗЛОВ MPIAIUoallv() >

КОНЕЦ СЧЕТА? (ДОСТИГНУТА НЕОБХОДИМАЯ \ ТОЧНОСТЬ)

НАЧАЛО

СЧИТЫВАНИЕ ФАЙЛОВ УЗЛОВ (.node) ТЕТРАЭДРОВ (.е\е)/

СБОР УЗЛОВ СЕТКИ НА КОРНЕВОМ ЯДРЕ MPIGatherv() /

ИТЕРАЦИЯ КГД АЛГОРИТМА

КОНЕЦ СЧЕТА? ^ (ДОСТИГНУТА НЕОБХОДИМАЯ ТОЧНОСТЬ}

КОРНЕВОЕ ЯДРО? MPICommrank()

СОХРАНЕНИЕ РЕЗУЛЬТАТОВ (.node, .ele) .

СОХРАНЕНИЕ РЕЗУЛЬТАТОВ (.node, ,ele) .

КОНЕЦ

КОНЕЦ

Рис. 3.1: Упрощенная блок-схема исходной последовательной программы моделирования течения газа (а) и ее параллельной версии (б) газодинамические величины) и набор тетраэдральных ячеек, определяющих связи между соседними вершинами. Подготовка расчетной сетки для параллельных расчетов заключается в разделении вычислительной работы между различными вычислительными узлами в соответствии с принципами геометрического параллелизма. То есть каждый вычислитель должен проводить расчеты газодинамических уравнений только на своем локальном множестве узлов сетки. Конечно-разностная схема расчета ((2.1)-(2.4) и (2.11)—(2.15)) такова, что для вычисления следующего по времени слоя в любом узле сетки требуются значения газодинамических величин во всех окружающих узлах. В связи с этим на каждой шаге по времени требуются обмены данными (средствами MPI) для передачи значений газодинамических величин в граничных узлах между каждой парой вычислителей, обрабатывающих соседние подобласти сетки. Представим расчетную сетку в виде графа, узлами которого являются узлы сетки, а ребрами - ребра тетраэдральных ячеек. Тогда общий объем передаваемой информации пропорционален количеству ребер, соединяющих вершины подграфов, соответствующих различным вычислителям. Затраты времени на передачу информации по сети средствами MPI (для самых распространенных сетевых решений - Ethernet, Myrinet, Infiniband) можно оценить следующим образом: в первом приближении они пропорциональны максимальному количеству ребер, исходящих из одного подграфа к узлам других подграфов. Следовательно, для подготовки расчетной сетки для параллельных расчетов требуется разбить исходный граф сетки на к подграфов таким образом, чтобы по возможности минимизировать время передачи данных (здесь к - количество процессорных ядер). Для статической балансировки вычислительной нагрузки необходимо, чтобы количество вершин в полученных подграфах было примерно одинаковым.

Расчетная сетка представляет собой два множества. Первое множество составляет набор узлов, в котором у каждого узла есть уникальный идентификатор и набор связанных с узлом переменных, соответствующих газодинамическим параметрам. Второе множество представляет собой набор тетраэдров, каждый из которых однозначно определяется четырьмя идентификаторами узлов. В таком виде расчетная сетка хранится в файлах и в оперативной памяти. Для представления ее в виде графа используется функция библиотеки METIS [58, 59] METISMeshToNodal(). Результатом работы этой функции является массив, в котором для каждого узла сетки перечислены все узлы, связанные с ним ребрами тетраэдральных ячеек.

Непосредственно разбиение графа осуществляется с помощью функции METISPartGraphKway(). Данная функция делит граф на к подграфов с равным количеством узлов, используя многоуровневый алгоритм [60, 61, 62, 63, 64, 65, 66, 67, 68, 69]. Разбиение производится таким образом, чтобы минимизировать общее количество граничных ребер. В процессе работы функция заполняет массив размерности N, где каждой вершине графа поставлено в соответствие число от 0 до (к — 1); здесь N - количество узлов сетки, к - количество процессоров, выделенных для расчета.

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

3.1.2 Построение схемы пересылки данных

В соответствии с построенными ранее численными схемами (2.1)-(2.4) и (2.11)—(2.15) для вычисления следующего временного слоя в любом узле требуются значения газодинамических переменных во всех соседних узлах сетки. Под соседними понимаются узлы, связанные с данным узлом ребром одной или более ячеек. Таким образом, на каждом шаге по времени все вычислители должны передавать значения в узлах своей подобласти тем своим соседям, которым это необходимо для расчета. В то же время необходимо принять от соседей значения газодинамических переменных в узлах, необходимых для расчета на следующем временном слое. Список, по которому происходит прием и передача узловых переменных, строится на каждом вычислителе один раз перед началом расчета. Для этого используется представление расчетной сетки в виде графа. Далее для каждого узла А локальной подобласти необходимо установить, имеет ли он связь с узлами другой подобласти. Предположим, что такой узел существует, и обозначим его В. Подобласть, которой он принадлежит, обозначим как j. Тогда в список принимаемых узлов от j-го вычислителя добавляется идентификатор узла В, а в список для передачи на j-ыЙ вычислитель добавляется идентификатор узла А. В обоих случаях используются глобальные идентификаторы узлов из исходной расчетной сетки (без разбиения). Алгоритм обхода узлов подобласти таков, что списки приема и передачи по построению оказываются отсортированными по идентификаторам узлов. Это обеспечивает взаимное совпадение списков приема и передачи на всех вычислителях.

Для хранения списков передачи используется контейнер vector стандартной библиотеки STL языка С++, элементами которого являются контейнеры set: include <vector> #include <set> typedef set<int> nodeSetT; typedef vector<nodeSetT> commVectorT;

Рассмотрим пример графа сетки (рис. 3.2). Построим векторы приема

Рис. 3.2: Граф сетки, разбиение по 4 ядрам и передачи узлов для ядра 1. Вектор передаваемых узлов:

0 [6 7 8 9]

1 []

2 [10 11 12]

3 []

Вектор принимаемых узлов:

0 [3 4 5]

1 []

2 [13 14]

3 []

3.1.3 Обмен данными (библиотека MPI)

Библиотека MPI [32, 33] предназначена для разработки программ, состоящих из независимых процессов, выполняющих свои наборы команд со своим набором данных (MIMD). То есть последовательности инструкций, исполняемые каждым процессом, не обязательно должны быть идентичными. Между собой процессы обмениваются посредством вызовов коммуникационных функций MPI.

MPI предоставляет протокол передачи данных, не зависящий от языка программирования. В нем реализованы как коммуникации "точка-точка", так и коллективные взаимодействия. Главные достоинства MPI - это высокая производительность, масштабируемость и портируемость. Высокая производительность достигается за счет того, что каждая конкретная реализация MPI специально оптимизирована под ту архитектуру, на которой она работает. Портируемость заключается в том, что реализации MPI существуют в настоящее время практически для любых параллельных архитектур. Большинство реализаций MPI состоят из набора функций (API), предназначенных для вызова из программ, написанных на языках С, С++, Fortran и других, поддерживающих вызовы данных библиотек.

В MPI включены интерфейсы для обеспечения виртуальных топологий, синхронизации, передачи данных между набором процессов (распределенных по отдельным вычислительным модулям). Общий принцип функционирования MPI не зависит от языка программирования, различаются только интерфейсы библиотеки. Основным понятием в MPI является процесс. Процесс ассоциируется с процессором (или ядром в многоядерных архитектурах), и обычно назначение одного процесса на каждый процессор является оптимальным, с точки зрения производительности, решением. Распределение MPI процессов по процессорам происходит при запуске MPI программы специальным агентом (обычно это mpirun или mpiexec).

Функции MPI включают операции взаимного приема/передачи "точка-точка", выбора декартовой логической топологии или топологии, задаваемой графом, обмена данными между парами процессов (прием/передача), объединения частичных результатов вычислений (операции сбора и редуцирования), синхронизации вычислительных узлов, получения системной информации (количество запущенных MPI процессов, идентификатор узла, на котором работает процесс, соседние процессы в логической топологии). Операции "точка-точка" подразделяются на синхронные и асинхронные, буфферизированные и блокирующие, предоставляя большую свободу программисту в выборе семантики.

В параллельном алгоритме моделирования течения газа на основе КГД уравнений обмены данными в течение одной итерации численной схемы происходят в два этапа. Сначала передаются данные внутренних узлов сетки на новом временном слое. Это требуется для аппроксимации граничных условий в соответствии с формулой (2.20). Затем передаются значения переменных в граничных узлах. В обоих случаях используется одна и та же функция, основанная на вызове MPI MPIAlltoallv(). Эта функция осуществляет разброс/сбор данных по принципу "все ко всем". Расчетная программа предварительно заполняет временные буферы значениями всех передаваемых переменных в соответствии с построенным ранее вектором передаваемых узлов. После выполнения вызова MPI расчетная программа распределяет полученные данные из соседних подобластей в соответствии с вектором принимаемых узлов.

Во время второго обмена данными кроме значений в граничных узлах передаются величина максимальной скорости звука max с (Pi) для коррекции шага по времени в соответствии с формулой 2.10 и величина текущей невязки. Оба значения вычисляются по узлам локальной подобласти, и после обмена каждый MPI процесс вычисляет глобальные максимальные значения скорости звука и невязки.

3.1.4 Параллельный ввод-вывод данных

Рассмотрим две важные особенности параллельного ввода-вывода в расчетной программе - чтение расчетной сетки (и сохраненных промежуточных данных) и запись результата в файлы. Несмотря на то, что при использовании библиотеки MPI допускается произвольный параллельный доступ на чтение к файлам из различных MPI процессов, функция чтения данных исходной последовательной программы требует модификации. Это обусловлено тем, что при использовании сетевых файловых систем (что является обычной ситуацией для многопользовательских высокопроизводительных вычислительных систем) чтение файлов данных при параллельном доступе фактически осуществляется последовательно. Поэтому время считывания одного и того же объема информации из файла группой из N вычислителей может вырасти почти в N раз.

В больших кластерных системах ("МВС-lOOk", "СКИФ МГУ") для поддержки файловых операций обычно используются более медленные сетевые решения (Ethernet), чем для передачи в среде MPI (Infiniband, Myrinet). В частности, в системе "МВС-lOOk" файловая система работает по технологии Gigabit Ethernet с пропускной способностью всего 100 Мбайт/с, при том что скорость MPI обменов примерно в 14 раз выше. Поэтому для сокращения времени считывания данных из файлов целесообразно производить это лишь в одном MPI процессе и передавать считанные данные в остальные процессы средствами MPI. Для этого в параллельной программе моделирования сверхзвуковых течений используется функция MPLBcast(). Данная функция позволяет за один вызов передать всем указанным процессам одинаковый блок данных. Проведенные эксперименты показали, что при таком подходе в практических расчетах время считывания данных увеличивается не более чем на 7-10% по сравнению с последовательной программой.

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

3.2 Описание оборудования

Расчеты модельных задач и оценка эффективности параллельных алгоритмов проводились на трех параллельных вычислительных системах. Задачи небольшого объема (задачи в двухмерной постановке, а также задачи, где число узлов расчетной сетки не превышало 100 тысяч) решались на кластере "IMM2" в Институте Математического Моделирования РАН. Задачи большего объема решались на кластерах "МВС-100k" Межведомственного Супер компьютерного Центра РАН и "СКИФ МГУ" в научно-исследовательском вычислительном центре МГУ имени М.В.Ломоносова.

Параметры вычислительных комплексов на момент проведения исследований приведены в таблицах 3.1, 3.2, 3.3.

Заключение с

Приведем основные результаты, полученные в рамках диссертационной работы:

1. Разработан численный алгоритм решения модельной газодинамической задачи на основе КГД уравнений на нерегулярных сетках большой размерности.

2. Разработана и реализована методика задания областей сложной формы с помощью сплайнов.

3. Разработан и реализован генератор сеток различного типа: треугольных, тетраэдральных и гибридных (сочетающих ячейки треугольной и прямоугольной формы).

4. Разработан параллельный алгоритм решения модельной газодинамической задачи.

5. Разработан комплекс параллельных программ для МВС различной архитектуры, ориентированный на решение задач внешнего обтекания тел сложной формы сверхзвуковым потоком вязкого газа на нерегулярных сетках большого объема.

6. Проведено тестирование комплекса на МВС различной архитектуры и мощности.

7. Проведены расчеты модельных двух- и трехмерных задач на сетках различного типа.

8. Проведено сравнение результатов моделирования с известными экспериментальными и расчетными дпанными других авторов.

К сказанному добавим следующее. Аппроксимации отдельных уравнений с частными производными на нерегулярных сетках были известны из литературы, однако их применение к КГД уравнениям в случае гибридных и тетраэдральных ячеек было проведено автором диссертации впервые. Методы аппроксимации границ тел с помощью сплайнов различного типа также были давно известны, однако в рамках данной работы потребовалась их специфическая реализация, которую автор выполнил самостоятельно. Разработанный генератор сеток в случае двумерных областей включает оригинальные программы, разработанные автором диссертации, а в случае трехмерных — использует в качестве базы пакет TetGen. Тестирование эффективности представленного параллельного алгоритма проводилось на вычислительных системах различного масштаба, и во всех случаях программный комплекс демонстрирует хороший уровень масштабируемости. Проведенные численные исследования подтвердили корректность используемых численных методов и применимость подхода для целого ряда практических задач.

1. В.А. Парфенов: Возвращение из космоса. Воениздат, 1961.

2. И.И. Ирбаева: Исследование границ начала и конца турбулентно-ламинарного перехода в пограничном слое в летном эксперименте при выведении ракет космического назначения. Ученые записки ЦАГИ, 36(1-2):70-78, 2005.

3. Самарский А.А., Михайлов А.П.: Математическое моделирование: Идеи. Методы. Примеры. Наука. Физматлит, Москва, 1997.

4. Gnoffo, P.A., K.J. Weilmuenster, Н.Н. Hamilton II, D.R. Olynick, and Е. Venkatapathy: Computational aerothermodynamic design issues for hypersonic vehicles. Journal of Spacecraft and Rockets, 36(l):21-43, 1999.

5. Елизарова Т.Г.: Квазигазодинамические уравнения и методы расчета вязких течений. Научный мир, Москва, 2007.

6. Шеретов Ю.В.: Анализ устойчивости модифицированной кинетически-согласованной схемы в акустическим приближении. В Применение функционального анализа в теории приближений, страницы 147-160. Тверской государственный университет, Тверь, 2004.

7. Елизарова Т.Г., Соколова М.Е.: Численный алгоритм расчета сверхзвуковых течений, основанный на квазигазодинамических уравнениях. Вестник Московского университета, серия 3. Физика. Астрономия, (1):10-15, 2004.

8. Кудряшова Т.А., Поляков С.В., Кононов Э.М.: Расчёт поля радиационного излучения газа вокруг спускаемого аппарата. Математическое моделирование, 20(10):63-74, 2008.

9. Graur I.A., Elizarova T.G., Lengrand J.С.: Investigation of shock configurations induced by a plume impinging upon a perpendicular plate. In Proc. of the 23rd Int. Symp. on Shock Waves, Fort Worth, Texas, US, 2001.

10. Елизарова Т.Г., Соколова М.Е.: Численный алгоритм расчета сверхзвуковых течений, основанный на квазигазодинамических уравнениях. Вестник Московского университета, серия 3. Физика. Астрономия, (1):10—15, 2004.

11. Елизарова Т.Г., Серегин В.В.: Аппроксимация квазигазодинамических уравнений на треугольных сетках. Вестник Московского университета, серия 3. Физика. Астрономия, (4): 15-18, 2005.

12. Четверушкин Б.Н.: Кинетические схемы и квазигазодинамическая система уравнений. Макс Пресс, Москва, 2004.

13. Граур И.А.: Алгоритмы численного решения квазигазодинамических уравнений. Журнал вычислительной математики и математической физики, 39(8): 1356-1371, 1999.

14. Chetverushkin B.N., Shilnikov E.V.: Kinetic finite-difference schemes and quasigasdynamic equation system as the physical model for gasdynamic flow description. In Proc. Ill Asian Comput. Fluid Dynamics Conf., pages 243-248, Bangalore, India, 1998.

15. Абалакин И.В., Жохова А.В.: Кинетически согласованные разностные схемы в обобщенной системе координат в двухмерном случае: Препринт №3. Диалог-МГУ, Москва, 1997.

16. Шильников Е.В., Шумков М.А.: 3D моделирование сверхзвукового течения вязкого газа вблизи каверны. Математическое моделирование, 11(10):17-30, 1999.

17. Шильников Е.В., Шумков М.А.: Моделирование трехмерных нестационарных течений газа на МВС с распределенной памятью. Математическое моделирование, 13(4):35-46, 2001.

18. Жохова А.В., Четверушкин Б.Н.: Моделирование нестационарных газодинамических течений. Математическое моделирование, 14(4):35-44, 2002.

19. Абалакин И.В., Жохова А.В., Четверушкин Б.Н.: Кинетически согласованные разностные схемы на нерегулярных сетках. Математическое моделирование, 9(7):44-53, 1997.

20. Абалакин И.В., Жохова А.В., Четверушкин Б.Н.: Кинетически согласованный алгоритм для расчета газодинамических течений на треугольных сетках. Математическое моделирование, 10(4):51-60, 1998.

21. А.В. Горобец: Масштабируемый алгоритм для моделирования несжимаемых течений на параллельных системах. Математическое моделирование, 19(10):105—128, 2007.

22. William L. Kleb, К. James Weilmuenster: A high angle of attack invis-cid shuttle orbiter computation. Technical Report NASA TM 107606, 1992.

23. Message passing interface (mpi) forum home page. http://www.mpi-forum.org/.

24. Корнеев В.Д.: Параллельное программирование в MPI. Институт компьютерных исследований, Москва-Ижевск, 2003.

25. The орептр api specification for parallel programming. http://www.openmp.org/.

26. Tetgen: A quality tetrahedral mesh generator. http://tetgen.berlios.de/.

27. Tecplot cfd post processing, plotting, graphing & data visualization software.http://www.tecplot.com/.

28. Попов Е.И.: Спускаемые аппараты. Новое в жизни, науке и технике. Космонавтика, астрономия, (4), 1985.38. 3d meshes research database by eric saltel inria gamma group. http://www-с.inria.fr/gamma/gamma.php.

29. Лойцянский JI.Г.: Механика жидкости и газа: Учеб. Для вузов -1-е изд., испр. Дрофа, Москва, 2003.

30. Dye, W. Н. and Т. Polek: Results of pressure distribution tests on a 0.010-scale space shuttle orbiter (61-0) in the nasa/arc 3.5-foot hypersonic wind tunnel (test oh38). Technical Report NASA CR-144584, 1975.

31. Fay, J. A. and F. R. Riddell: Theory of stagnation point heat transfer in dissociated air. Journal of Aeronautical Science, 25(l):73-85, 1958.

32. Ames Research Center Staff: Equations, tables, and charts for compressible flow. Technical Report NACA Report 1135, 1953.

33. Елизарова Т.Г., Хохлов А.А.: Квазигазодинамические уравнения для течений газа с внешними источниками тепла. Вестник Московского университета, серия 3. Физика. Астрономия, (2):22-25, 2007.

34. Елизарова Т.Г., Четверушкин Б.Н.: Кинетический алгоритм для расчета газодинамических течений. Журнал вычислительной математики и математической физики, 25(10):1526—1533, 1985.

35. Шеретов Ю.В.: Математическое моделирование течений жидкости и газа на основе квазигидродинамических и квазигазодинамических уравнений. Тверской государственный университет, Тверь, 2000.

36. Елизарова Т.Г., Шеретов Ю.В.: Теоретическое и численное исследование квазигазодинамических и квазигидродинамических уравнений. Журнал вычислительной математики и математической физики, 41(2):239-255, 2001.

37. Елизарова Т.Г., Соколова М.Е., Шеретов Ю.В.: Квазигазодинамические уравнения и численное моделирование течений вязкого газа. Журнал вычислительной математики и математической физики, 45(3):544—555, 2005.

38. Ильгамов М.А., Гильманов А.Н.: Неотражающие условия на границах расчетной области. Физматлит, Москва, 2003.

39. А.А. Самарский: Введение в численные методы. Наука, Москва, 1987.

40. Т. Кормен, Ч. Лейзерсон, Р. Ривест: Алгоритмы: построение и анализ. МЦНМО, Москва, 2001.

41. Domain decomposition methods. http://www.ddm.org/.

42. I.V. Abalakin, A. Dervieux, Т.К. Kozubskaya: High accuracy finite volume method for solving nonlinear aeroacoustics problems on unstructured meshes. Chinese Journal of Aeroanautics, 19(2), 2006.

43. H. Edelsbrunner, N.R. Shah: Incremental topological flipping works for regular triangulations. Algorithmica, 15:223-241, 1996.

44. H. Si, K. Gaertner: Meshing piecewise linear complexes by constrained delaunay tetrahedralizations. In Proceedings of the 14th International Meshing Roundtable, pages 147-163, 2005.

45. Si, H.: On refinement of constrained delaunay tetrahedralizations. In Proceedings of the 15th International Meshing Roundtable, 2006.

46. Карамзин Ю.Н., Попов И.В., Поляков С.В.: Разностные методы решения задач механики сплошной среды на неструктурированных треугольных и тетраэдральных сетках. Математическое моделирование, 15(11), 2003.

47. Shewchuk, Jonathan Richard: Delaunay Refinement Mesh Generation. PhD thesis, School of Computer Science, Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1997. Technical Report CMU-CS-97-137http://www.cs.emu.edu/~jrs/jrspapers.html.

48. Metis family of multilevel partitioning algorithms. http://glaros.dtc.umn.edu/gkhome/views/metis.

49. G. Karypis, K. Schloegel, V. Kumar: Parmetis: Parallel graph partitioning and sparse matrix ordering library. Technical report, Department of Computer Science, University of Minnesota, Minneapolis, MN, 1998.

50. B. Hendrickson, R. Leland: A multilevel algorithm for partitioning graphs. In Supercomputing '95 Proceedings, San Diego, CA, 1995.

51. B. Hendrickson, R. Leland: An improved spectral graph partitioning algorithm for mapping parallel computations. SIAM Journal on Scientific Computing, 16(2):452-469, 1995. *

52. B. Hendrickson, R. Leland: An improved spectral load balancing method. In Proceedings of the 6th SIAM Conference on Parallel Processing for Scientific Computing, pages 953-961. Society for Industrial and Applied Mathematics, 1993.

53. B. Hendrickson, R. Leland: Multidimensional spectral load balancing. Technical report SAND 93-0074, Sandia National Laboratories, Albuquerque, NM, 1993.

54. G. Karypis, V. Kumar: Multilevel algorithms for multi-containt graph partitioning. Technical report TR 98-019, Department of Computer Science, University of Minnesota, Minneapolis, MN, 1998.

55. G. Karypis, V. Kumar: Multilevel k-way partitioning scheme for irregular graphs. Journal of Parallel and Distributed Computing, 48(1) :96-129, 1998.

56. G. Karypis, V. Kumar: Parallel multilevel k-way partitioning scheme for irregular graphs. Siam Review, 41(2):278-300, 1999.

57. G. Karypis, V. Kumar: A fast and high quality multilevel scheme for partitioning irregular graphs. Technical report CORR 95-035, Department of Computer Science, University of Minnesota, Minneapolis, MN, 1995.

58. G. Karypis, V. Kumar: Analysis of multilevel graph partitioning. In Proceedings of the 7th Supercomputing Conference, 1995.

59. Amdahl, Gene: Validity of the single processor approach to achieving large-scale computing capabilities. In AFIPS Conference Proceedings, number 30, pages 483-485, 1967.

60. Boost С++ libraries. http://www.boost.org/.

61. Prey, Pascal J: Medit : An interactive mesh visualization software. Technical report RT-0253, INRIA-Rocquencourt, 2001.

62. Gnoffo, P. A.: An upwind-biased, point-implicit relaxation algorithm for viscous, compressible perfect-gas flows. Technical Report NASA TP 2953, 1990.

63. K. J.Weilmuenster and H. H. Hamilton, II: Calculations of inviscid flowover shuttle-like vehicles at high angles of attack and comparisons with experimental data. Technical Report TP-2103, NASA, 1983.1. U3