CC BY
93
УДК 519.6
УЧЕБНО-ЭКСПЕРИМЕНТАЛЬНЫЙ ВЫЧИСЛИТЕЛЬНЫЙ КЛАСТЕР. Ч. 2. ПРИМЕРЫ РЕШЕНИЯ ЗАДАЧ
В.В. Лукин, И.К. Марчевский, В.С. Морева, А.Ю. Попов, К.Л. Шаповалов, Г.А. Щеглов
МГТУ им. Н.Э. Баумана, Москва
e-mail: wlukin@gmail.com; iliamarchevsky@mail.ru;
morevavs@rambler.ru; georg@energomen.ru
Представлены разработанные параллельные численные алгоритмы решения задач магнитной гидродинамики методом RKDG и некоторых классов задач аэрогидродинамики методами вихревых элементов. Приведены примеры решения задач на Учебно-экспериментальном вычислительном кластере кафедры "Прикладная математика" МГТУ им. Н.Э.Баумана. Исследована эффективность распараллеливания вычислений в пакете OpenFOAM и сформулированы рекомендации по повышению эффективности работы кластера. Рассмотрен вопрос разработки учебного курса по параллельньш вычислениям, включенного в учебный план кафедры.
Ключевые слова: параллельные вычисления, вычислительный кластер, математическое моделирование, магнитная гидродинамика, метод вихревых элементов, пакет OpenFOAM.
COMPUTING CLUSTER FOR TRAINING
AND EXPERIMENTS. PART 2. EXAMPLES OF SOLVING PROBLEMS
V.V. Lukin, I.K. Marchevskii, V.S. Moreva, A.Yu. Popov, K.L. Shapovalov, G.A. Shcheglov
Bauman Moscow State Technical University, Moscow e-mail: vvlukin@gmail.com; iliamarchevsky@mail.ru morevavs@rambler.ru; georg@energomen.ru
The developed parallel numerical algorithms are presented for solving problems of magnetic hydrodynamics by the RKDG method and certain classes of problems of aerodynamics by vortex element methods. Examples of solving the problems at the computing cluster for training and experiments of the "Applied Mathematics" department of the Bauman Moscow State Technical University are given. The efficiency of parallelizing computations in the OpenFOAM software package is studied, and the recommendations on improvement of the cluster operation efficiency are formulated. An issue of development of the learning course on parallel computing included in the department curriculum is considered.
Keywords: parallel computing, computing cluster, mathematical modeling, magnetic hydrodynamics, vortex element method, OpenFOAM software package.
На кафедре ФН-2 "Прикладная математика" МГТУ им. Н.Э. Баумана в 2011г. с целью проведения исследований в области параллельных алгоритмов решения задач математической физики создан Учебно-экспериментальный вычислительный кластер. В работе [1] приведены основные технические характеристики созданной высокопроизводительной системы, причем при описании использована специальная терминология, позволившая корректно разделить логический и аппаратный подходы к описанию кластера.
Кратко напомним основные понятия введенной терминологии:
• кластер — это объединение логической и аппаратной структуры вычислительной системы;
• вычислительное ядро — единица логической структуры кластера, элементарный вычислитель, способный исполнять последовательную программу (вычислительное ядро следует отличать от процессорного ядра, которое можно рассматривать как элемент аппаратной структуры);
• вычислительный модуль — элемент логической структуры кластера, представляющий собой одно или несколько вычислительных ядер, объединенных общим полем памяти;
• узел — элемент аппаратной структуры кластера, объединяющий один или несколько модулей в сервер.
Учебно-экспериментальный вычислительный кластер состоит из одного управляющего и восьми вычислительных узлов на базе 4-ядерных процессоров Intel Core i7-920, объединенных сетью Gigabit Ethernet, общий объем оперативной памяти — 100 Гб. Максимальная производительность кластера, измеренная по тесту Linpack Benchmark, 262,2 Гфлопс [1], что составляет 79% пиковой производительности.
Основные цели создания кластера — это решение ресурсоемких прикладных задач математической физики в рамках научных направлений, развиваемых на кафедре ФН-2 "Прикладная математика" МГТУ им. Н.Э. Баумана, а также обучение студентов и аспирантов подходам и приемам разработки эффективных параллельных вычислительных алгоритмов. В настоящей статье расмотрены некоторые алгоритмы для решения задач аэрогидродинамики, магнитной гидродинамики и тепломассопереноса. В 2012 г. авторами был организован и проведен учебно-научный семинар "Параллельные вычислительные технологии" для студентов и аспирантов МГТУ им. Н.Э. Баумана, на котором рассмотрены базовые вопросы организации высокопроизводительных вычислений, архитектуры многопроцессорных вычислительных комплексов и основные технологии параллельного программирования для различных классов вычислительных систем, такие как OpenMP, MPI, nVidia CUDA. Описание курса приведено в настоящей работе.
Решение прикладных задач. Ниже рассмотрено несколько параллельных алгоритмов решения практических задач, описываемых математическими моделями магнитной гидродинамики, динамики идеальной и вязкой несжимаемой среды, тепломассопереноса в несжимаемой среде. Для каждой задачи приведена математическая постановка, кратко описан используемый численный метод и способ его параллельной
реализации. Для оценки эффективности распараллеливания на кластере проведены тестовые расчеты. За основной критерий эффективности параллельного алгоритма принято получаемое ускорение вычислений при вовлечении в расчет различного числа вычислительных модулей. При этом под ускорением понимается отношение времени решения задачи с помощью последовательного алгоритма ко времени ее решения параллельным алгоритмом [2].
Параллельная реализация метода RKDG для решения задач идеальной магнитной гидродинамики. Уравнения идеальной магнитной гидродинамики (МГД) описывают большое число процессов, происходящих как в технических системах, так и, например, в астрофизических объектах. Особый интерес представляет процесс развития неустойчивости течения в окрестности поверхности разрыва, для численного моделирования которого целесообразно использовать методы второго и более высоких порядков точности. Один из таких методов — RKDG (Runge-Kutta Discontinuous Galerkin), поскольку он обеспечивает высокое разрешение разрывов и позволяет повышать порядок точности на неструктурированных сетках без увеличения шаблона аппроксимации [3]. Авторами разработан и апробирован программный комплекс PLASMA, позволяющий решать двумерные уравнения идеальной МГД с помощью метода RKDG.
Математическая модель рассматриваемого класса задач содержит уравнения, выражающие законы сохранения массы, импульса и энергии, а также уравнение Фарадея для магнитной индукции [4]:
Здесь р — плотность газа, р — давление, е — полная энергия единицы объема, V = (у\, г>2, ^з)т — вектор скорости, В = (В\, В2, £>з)т -
О
(1)
где
вектор магнитной индукции. Из закона Фарадея следует условие без-дивергентности магнитного поля сЦу В = 0.
Система (1) замыкается уравнением состояния совершенного газа р = (у — 1)ре, где 7 — показатель адиабаты, е — удельная внутренняя
энергия. Отсюда е =
Р
/ЛГ
-. Можно показать, что эта
7-1 2 8тг система гиперболическая [4] и, следовательно, существует разложение
Ж
ди
(2)
где Ьг и Лг — ортогональные матрицы, составленные из левых и правых собственных векторов, Л^ — действительная диагональная матрица, составленная из собственных чисел матрицы А^
Для численного решения системы (1) на треугольной неструктурированной сетке применяется разрывный метод Галеркина [3], в соответствии с которым решение и представляется в виде
где Л^ — число ячеек сетки; Щ — зависящее от порядка аппроксимации число базисных функций, определенных в ячейке; ъ^г)}^! —
ортонормированная система базисных функций (полиномов), определенных в 7'-й ячейке; у,/, (£) — вектор коэффициентов при базисных функциях.
После домножения обеих частей (1) на систему пробных функ-
интегрирования по ячейкам и некоторых пре-
ций У {^А
образований приходим к системе обыкновенных дифференциальных уравнений
dy.
dt
з
§
F* (Vi, Ути,
dS-
Sa
ъ
Fo
где а — номер ребра ячейки, а) — номер соседней с ней ячейки, — численный поток, зависящий от решений в
ячейках, примыкающих к ребру а, и определяемый из решения соответствующей задачи Римана о распаде разрыва.
Система (З) дополняется начальным условием вида
где
начальное распределение консервативных перемен-
ных.
Для приближенного решения задачи Римана для системы МГД уравнений используется ряд методов с выделением разрывов, такие как НЬЬ, НЬЬС, НЬЬБ [4, 5]. Метод НЬЬБ позволяет получить наилучшее разрешение разрывов, но во многих случаях приводит к возникновению численной немонотонности в решении. В этих случаях используется более диссипативный численный поток НЬЬ. При расчете интегралы заменяются квадратурными формулами, точность которых согласована с порядком метода. Решение задачи Коши (3), (4) проводится численным интегрированием явным методом Рунге-Кутты, шаг по времени Д£ определяется динамически из условия устойчи-
АЬ
вости Куранта-Фридрихса-Леви: тах|Аг|-— < 1, где /¿,мт — длина
наименьшего ребра ячейки в сетке, A¿ — собственные числа системы (1) на предыдущем шаге по времени.
Для сохранения монотонности решения, а также для предотвращения появления в решении отрицательных значений плотности и давления в случае применения метода второго порядка точности используются лимитеры — ограничители наклона решения в ячейке [4]. В программе PLASMA реализован "TVB minmod" лимитер для кусочно-линейных функций [3, 4]. Так как процедура лимитирования магнитного поля приводит к накоплению численного магнитного заряда, она применяется только к величинам, используемым при вычислении потоков на границах ячеек. Лимитирование проводится для локальных характеристических переменных.
Для поддержания бездивергентности магнитного поля происходит пересчет магнитных составляющих потоков, получаемых в результате решения задачи Римана [6]. Аппроксимация граничных условий проведена методом фиктивной ячейки.
Нахождение решения системы (1) методом RKDG — вычислительно сложная задача и при достаточно больших размерах сетки требует применения параллельных вычислительных технологий. Параллельный вариант программы реализован при помощи технологии MPI (Message Passing Interface) [7] и основан на разбиении исходной области на подобласти, решение в каждой из которых рассчитывается отдельным вычислительным модулем. В этом случае алгоритм решения
min
Рис. 1. Пример разбиения сложной области на 16 подобластей при помощи программы METIS
состоит из параллельного нахождения решения на текущем временном слое в заданных подобластях и последующего, реализованного с помощью средств MPI, обмена "контактными" граничными данными между модулями, обрабатывающими соседние подобласти.
Для достижения максимальной эффективности распараллеливания вычислений необходимо минимизировать количество граничных данных, пересылаемых между вычислительными модулями (число "контактных" ячеек в каждой подобласти), и сбалансировать число расчетных ячеек в каждой подобласти. Для приближенного нахождения подобного разбиения использована программа METIS (рис. 1).
На основании полученного разбиения главный (нулевой) вычислительный модуль определяет и рассылает массивы, содержащие локальную нумерацию ячеек в подобластях и номера "контактных" ячеек. Каждый вычислительный модуль хранит в памяти информацию только о своей подобласти и ее границе, а также о расположении соседних подобластей. На основании разосланных массивов каждый из модулей проводит инициализацию отложенного приема и передачи "контактных" данных. Далее вычислительные модули проводят расчет решения независимо, выполняя на каждом шаге по времени обмен граничными данными и определение величины нового шага по времени на основании максимальной по всем подобластям скорости распространения возмущений. Периодически главный вычислительный модуль проводит сбор и сохранение решения на текущем временном слое.
В качестве тестовой рассмотрена задача о вращении цилиндра в покоящейся среде [8, 9]. Расчетная область представляет собой единичный квадрат. В центре области находится вращающийся плотный диск радиусом 0,1, который постепенно закручивает окружающее его замагниченное вещество. Начальные условия задачи:
(1, О, О, 0, 1, 5%]%, г'>'о,115'
0,5)=
На границе области задано условие свободного вытекания вещества. Показатель адиабаты 7 = 5/3, что соответствует идеальному двухатомному газу, число Куранта С = 0,4. Расчеты проведены на треугольных сетках, содержащих 150 000 и 550 000 ячеек. На рис.2 представлены распределения модуля вектора скорости в момент времени t = 0,15, полученные разрывными методами Галеркина первого и второго порядка точности с использованием численного потока HLLD (для грубой сетки). Решение методом второго порядка приводит к меньшему размытию разрывов как на внешней границе возмущенного вещества, так и в области повышенной плотности.
На рис. 3 представлен график ускорения параллельной программы при решении указанной задачи на учебно-экспериментальном вычислительном кластере, а также графики линейного (идеального) ускорения и максимального ускорения согласно закону Амдала.
Доля параллельного кода в описанном алгоритме для подробной сетки составляет 99,9%, что согласно закону Амдала ограничивает ускорение величиной 31,2 на 32 ядрах. Решение задачи на одном ядре занимает 234,3 мин, на 32 ядрах — 10,1 мин (ускорение в 23,2 раза). Существенное расхождение реального ускорения с оценкой по закону Амдала объясняется главным образом не влиянием межмодульных
Рис. 2. Распределения модуля вектора скорости в задаче о вращении цилиндра, полученные методами ККОС первого (я) и второго (в) порядков
Ускорение
- Линейное ускорение --Ускорение по Амдалу
/V"
—•— Реальное ускорение (подробная сетка)
— Реальное ускорение (грубая сетка)
Число модулей
4 8 12 16 20 24 28 32
Рис.3. Графики ускорения расчетов в задаче о вращении цилиндра в покоящейся среде
обменов, а спецификой аппаратной структуры кластера: вычислительные модули сгруппированы по четыре в рамках восьми вычислительных узлов, что приводит к конкуренции запросов модулей одного узла к оперативной памяти, а следовательно — к большому числу кэш-промахов. Это особенно критично для текущей реализации метода решения уравнений МГД второго порядка точности, поскольку каждой ячейке сетки соответствует как минимум 32 числа с плавающей точкой. Расчеты показали, что процедура вычисления потоков через ребра на четырех ядрах ускоряется в 3,8 раза, в то время как копирование решения с предыдущего слоя на новый — только в 2,1 раза.
Параллельная реализация метода вихревых элементов. Метод вихревых элементов (МВЭ) является весьма эффективным численным методом моделирования двумерных течений несжимаемой среды, и в настоящее время он интенсивно развивается как в нашей стране, так и за рубежом [10, 11]. Распределение завихренности в МВЭ моделируется набором вихревых элементов (ВЭ), для каждого из которых задается положение в пространстве и интенсивность. Скорость среды в любой точке течения вычисляется по закону Био-Савара по известным характеристикам ВЭ, а давление определяется при помощи аналога интеграла Коши-Лагранжа.
Возможности метода достаточно широки: он, с одной стороны, позволяет исследовать фундаментальные вопросы гидромеханики вплоть до эффектов, которые невозможно или чрезвычайно трудно смоделировать при использовании других численных методов [11], а с другой стороны — дает возможность решать многие инженерные задачи с приемлемой для практических приложений точностью при достаточно низких затратах машинного времени по сравнению с сеточными методами. Одна из актуальных задач — определение зависимостей аэродинамических коэффициентов профиля Сха, Суа и Ст — коэффици-
ентов лобового сопротивления, подъемной силы и аэродинамического момента соответственно — от угла атаки.
Для получения стационарных значений Сха, С,„, и Ст для фиксированного угла атаки профиля требуется решить нестационарную задачу о моделировании обтекания неподвижного профиля, установленного под заданным углом, а затем осреднить вычисленные нестационарные значения коэффициентов по большому промежутку времени [12]. Для получения достаточно точных зависимостей аэродинамических характеристик от угла атаки может потребоваться провести много расчетов. Скорость решения всей задачи, т.е. построения зависимостей Сха{а), Суа{р) и Ст(а) для определенного диапазона углов атаки, становится особенно важной, когда требуется провести анализ большого числа возможных вариантов, например при решении задачи оптимизации формы профиля. Аналогичная ситуация возникает при проведении вычислительных экспериментов, чрезвычайно важных при разработке новых или модификации существующих расчетных схем МВЭ. К этому же классу относятся "методические" задачи подбора оптимальных параметров расчетной схемы численного метода.
Течение несжимаемой среды постоянной плотности р описывается уравнениями неразрывности и Навье-Стокса:
дУ _ [р\
(V-
дЬ ' ■ \р/ ■
где V — скорость среды, р — давление, и — коэффициент кинематической вязкости. На бесконечности выполняется граничное условие затухания возмущений
а на поверхности неподвижного обтекаемого профиля К — условие прилипания
Особенно эффективен МВЭ именно при решении внешних задач обтекания, поскольку на каждом шаге расчета по времени выполнение уравнения неразрывности и граничного условия на бесконечности происходит автоматически. Граничное условие прилипания обеспечивается генерацией ВЭ на поверхности обтекаемых профилей, а течение среды моделируется движением имеющихся ВЭ. Интенсивности рождаемых ВЭ находятся из решения соответствующей системы линейных алгебраических уравнений, аппроксимирующей граничное условие на профиле, а движение ВЭ описывается системой обыкновенных дифференциальных уравнений, правые части которых представляют
Л*."
Рис. 4. Обтекание кругового цилиндра (число Рейнольдса Яс = 1000)
собой скорости движения ВЭ. На рис. 4 приведен пример расчета обтекания кругового цилиндра с образованием за ним вихревой дорожки Кармана.
Процедура определения скоростей вихревых элементов является наиболее затратной с вычислительной точки зрения, особенно при увеличении числа ВЭ, поскольку для этого приходится учитывать взаимное влияние всех пар ВЭ. Несколько менее трудоемкими, но тем не менее достаточно затратными являются операции реструктуризации вихревой пелены и вычисления аэродинамических нагрузок, действующих на профиль.
Таким образом, решение серии задач расчета обтекания профилей под различными углами атаки требует выполнения большого объема вычислительной работы, при этом следует учитывать однотипность и независимость отдельных задач. В каждой из задач с помощью алгоритма МВЭ моделируется обтекание профиля под фиксированным углом атаки и определяются аэродинамические нагрузки, действующие на него. В силу независимости этих задач наиболее очевидным является внешнее распараллеливание — одновременное решение нескольких задач. Использование большого числа вычислительных модулей дает возможность проводить несколько расчетов одновременно; по мере завершения отдельных расчетов и освобождения соответствующих модулей на них запускается процесс решения следующих задач.
В разработанном авторами программном комплекс РОЬАКА реализовано также внутреннее распараллеливание — параллельное выполнение наиболее трудоемких этапов алгоритма внутри решения одной задачи расчета обтекания профиля. Использование данного механизма позволяет оптимизировать загрузку вычислительных модулей при решении большой серии задач, а также сократить временные затраты при решении малой серии задач или единственной задачи, которая может иметь высокую вычислительную сложность. Возможно использование обоих типов параллелизма одновременно, причем число вычислительных модулей, выделяемых под решение различных
задач расчета обтекания профиля под конкретными углами атаки, может выбираться пользователем произвольно. Это связано с тем, что для одного и того же профиля вычислительная сложность моделирования его обтекания может существенно зависеть от угла атаки.
Наиболее затратными с вычислительной точки зрения являются следующие операции алгоритма МВЭ.
Операция 1. Вычисление скоростей ВЭ как результата их взаимного влияния друг на друга.
Операция 2. Реструктуризация вихревой пелены, приводящая к уменьшению числа ВЭ за счет объединения близкорасположенных ВЭ или исключения из расчетной схемы ВЭ, удалившихся от обтекаемого профиля на значительное расстояние.
Операция 3. Вычисление нагрузок, которые действуют на профиль со стороны потока, по рассчитанному полю скоростей и известному распределению завихренности.
Проведены исследования эффективности использования параллельных алгоритмов при решении плоских задач МВЭ (внутреннего распараллеливания). На рис. 5 представлены оценки вычислительной трудоемкости различных этапов последовательного алгоритма МВЭ. Наиболее трудоемким является вычисление скорости ВЭ. Другие две операции — реструктуризация вихревой пелены и вычисление аэродинамических нагрузок менее трудоемки, в то время как затратами машинного времени на прочие операции (например, решение системы линейных уравнений) можно пренебречь.
С помощью закона Амдала можно оценить максимальное ускорение, которое теоретически достижимо при распараллеливании только операции 1, операций 1, 2 или операций 1, 2, 3 алгоритма, например, для вычислительных систем с 4, 16 и 64 вычислительными модулями. Результаты представлены в табл. 1.
Рис. 5. Трудоемкость операций в методе вихревых элементов
Таблица 1
Максимальное ускорение параллельного кода
Распараллеленные Доля параллельного Максимальное ускорение
операции кода 4 модуля 16 модулей 64 модуля
1 95,3 % 3,5 9,4 16,2
1, 2 98,2 % 3,8 12,6 30,0
1, 2, 3 99,5 % 3,9 14,9 48,7
Для вычислений с использованием небольшого числа вычислительных модулей максимальные значения ускорений близки. Но при использовании системы с десятками вычислительных модулей распараллеливание операций 2 и 3 дает значительный вклад в ускорение, несмотря на то что трудоемкость этих операций составляют вместе лишь малую долю трудоемкости операции 1 (приблизительно 4,4 %).
На рис. 6 для программного комплекса РОЬАКА представлена зависимость ускорения от числа вычислительных ядер при распараллеливании операций 1, 2 и 3 алгоритма, полученная при решении одной весьма сложной с вычислительной точки зрения задачи расчета обтекания профиля под фиксированным углом атаки. Для сравнения показано максимальное значение ускорения для случаев распараллеливания только операции 1, операций 1, 2 и операций 1, 2, 3.
Реальное ускорение при распараллеливании операций 1, 2, 3 близко к оценке по закону Амдала для распараллеливания только операций 1, 2 (сказывается наличие пересылок данных, затрат на синхронизацию вычислений и прочих "накладных расходов"), поэтому можно сделать вывод о том, что распараллеливание операции 3 важно для достижения высокой эффективности, хотя ее трудоемкость крайне мала — всего около 1,3%.
Ускорение
Распараллеливание операции 1
-о- Распараллеливание операций 1, 2
- ■ - Распараллеливание операций 1, 2, 3
—и— Реальное ускорение
О 8 16 24 32 40 48 56 64 72 80
Рис. 6. Реальное ускорение и оценки по закону Амдала (внутреннее распараллеливание)
Если требуется решить не одну задачу, пусть и весьма трудоемкую, а серию задач расчета обтекания профиля под различными углами атаки, эффективность применения параллельных алгоритмов (внешнего распараллеливания) может быть существенно выше. В качестве тестового примера рассматривалась серия из 91 задачи по моделированию обтекания крылового профиля под разными углами атаки; в каждом расчете выполнялось 5000 временных шагов. На рис. 7 показаны затраты машинного времени при решении этой серии задач с помощью программного комплекса POLARA на 1, 32 и 64 вычислительных модулях, а также зависимость ускорения вычислений от числа используемых ядер при проведении расчетов на вычислительных комплексах МВС-ЮОК (МСЦ РАН), СКИФ-МГУ "Чебышев" (МГУ им.М.В.Ломоносова) и Учебно-экспериментальном вычислительном кластере.
На всех вычислительных комплексах ускорение растет практически линейно с ростом числа задействованных модулей; при использовании 64 вычислительных модулей удается получить более, чем 50-кратное ускорение расчетов.
В рассмотренном выше примере для распараллеливания вычислений использовалась библиотека MPI как для внешнего, так и для внутреннего распараллеливания. Целесообразность ее применения для внешнего распараллеливания, т.е. для одновременного решения нескольких полностью независимых задач, не вызывает сомнений ввиду практически отсутствующего межмодульного обмена данными. В то же время при распараллеливании решения каждой задачи (внутрен-
Рис. 7. Ускорение при решении серии задач (внешнее распараллеливание)
нем) интенсивность обменов данными между задействованными в расчете вычислительными модулями становится существенной, поэтому для повышения эффективности работы программного комплекса была исследована возможность применения для этих целей технологии ОрепМР. С ее помощью были распараллелены те же три наиболее трудозатратные операции, что и при использовании МР1.
При проведении расчетов на вычислительных модулях с общей памятью на 2-8 вычислительных ядрах установлено, что выигрыш по времени, достигаемый за счет использования ОрепМР, не превышает нескольких процентов по сравнению с применением МР1. В связи с этим данная возможность практически не использовалась, поскольку применение MPI позволяет сделать программу универсальной — она может работать на системах как с общей, так и с распределенной памятью. При этом для целей внутреннего распараллеливания может использоваться произвольное число вычислительных модулей, которые могут находиться на различных узлах кластера. В последнем случае эффективность распараллеливания снизится, но такая возможность необходима для решения сложных задач, когда в расчетной схеме присутствуют десятки-сотни тысяч ВЭ.
Параллельные алгоритмы моделирования динамики вихревых структур. Метод вихревых элементов может быть обобщен на случай моделирования пространственных течений идеальной несжимаемой среды. Однако если в двумерном случае в качестве ВЭ всегда используется изолированный круглый вихрь (вихрь Рэнкина), то для решения пространственных задач разработано достаточно много моделей различных ВЭ, каждая из которых имеет определенные достоинства и недостатки. Также существенную сложность представляет "сопряжение" вихревого слоя на поверхности обтекаемого тела с полем завихренности в потоке.
В работе [13] предложена модель симметричного вортона-отрезка, позволяющая с достаточной для практических целей точностью и приемлемыми затратами машинного времени решать широкий класс задач расчета внешнего обтекания тел и особенно эффективная при моделировании динамики вихревых структур. Она, в частности, позволяет моделировать такие эффекты как изгиб, удлинение и перезамыкание вихревых нитей. При решении этих задач область течения считается безграничной, в начальный момент времени в нее помещается вихревая структура, дальнейшее движение которой моделируется в процессе расчета. Данная вихревая структура моделируется множеством из N вортонов-отрезков, каждый из которых характеризуется двумя векторными параметрами — радиус-вектором вортона г, и вихревым вектором hi, а также скалярным параметром — интенсивностью вихревого отрезка Г\.
В процессе эволюции вихревой структуры интенсивности ВЭ остаются неизменными (I\ = const), а изменение векторов г, и h описывается системой обыкновенных дифференциальных уравнений
где Vi — скорость, индуцированная всеми вортонами-отрезками в месте расположения г-го вортона; [В,] — тензор деформации г-го вортона-отрезка, равный градиенту поля скоростей в месте расположения г-го вортона, горизонтальная черта означает усреднение соответствующей величины вдоль вортона-отрезка.
Для вычисления скоростей, индуцируемых вортонами-отрезками, а также их градиентов в работе [13] получены аналитические формулы, в то время как их интегрирование при вычислении средних значений проводится численно с использованием квадратурных формул Гаусса с тремя или пятью точками.
Для оценки точности моделирования эффектов изгиба и перезамыкания вихревых колец, составленных из симметричных вортонов-отрезков, рассмотрена задача об эволюции овального вихревого кольца с удлинением 7,5. На рис. 8 показаны результаты расчета эволюции кольца единичной циркуляции. Как видно из рисунка, при t ~ 2,0 происходит первое перезамыкание вихревых нитей с образованием двух колец из одного, а затем при t ~ 5,5 происходит второе перезамыкание нитей и два кольца объединяются в одно. Таким образом, использование симметричного вортона-отрезка позволяет устойчиво моделировать перезамыкание вихревых нитей с малой погрешностью.
На рис. 9 приведены графики ускорения расчета в зависимости от числа задействованных вортонов.
Эффективность распараллеливания при большом числе используемых вычислительных модулей существенно зависит от трудоемкости задачи. При наличии в расчетной схеме 20 ООО вортонов и более ускорение близко к линейному.
Параллельные вычисления в пакете OpenFOAM. Открытый свободно распространяемый пакет OpenFOAM применяется для решения широкого класса задач механики сплошной среды. Он представляет собой совокупность большого числа расчетных подпрограмм для решения различных классов уравнений в частных производных, а также специальных утилит для подготовки исходных данных, сопровождения вычислений и обработки результатов вычислений.
Отличительной особенностью OpenFOAM является его ориентация на выполнение вычислений на высокопроизводительных кластерных системах. Все расчетные подпрограммы и большая часть утилит
Рис. 8. Фазы процесса эволюции овального вихревого кольца
Ускорение
-*- 2000 вортонов -«-3500 вортонов -»«-30000 вортонов — Линейное ускорение
и-" Число модулей
4 8 12 16 20 24 28 32
Рис. 9. Графики ускорения расчетов при различном числе вортонов
пре- и постпроцессинга могут быть запущены в параллельном режиме с использованием технологии MPI. При этом процедура запуска параллельной версии мало отличается от последовательной; от пользователя требуется лишь указать число вычислительных модулей, которые будут задействованы в расчете, и способ разделения расчетной сетки на подобласти. Такое разбиение может задаваться как явным указанием границ подобластей, так и в автоматизированном режиме с использованием алгоритмов SCOTCH или METIS.
Отметим, что алгоритмы OpenFOAM не предполагают распараллеливания процедуры решения систем линейных алгебраических уравнений (СЛАУ), возникающих на каждом шаге расчета по времени. Это связано со сложностями разработки эффективных параллельных
Рис. 10. График ускорения расчетов при решении задач на сетке с 55 ООО ячеек
реализации современных итерационных численных методов, реализованных в данном пакете, — проекционных методов крыловского типа и многосеточных методов. На каждом вычислительном модуле для сво-еИ подобласти формируется отдельная СЛАУ, а согласование решения на межмодульных границах осуществляется итерационно.
Эффективность распараллеливания решения задачи в ОрепБОАМ зависит от многих факторов, прежде всего, от вычислительнои сложности самои задачи, под которои понимается время, затрачиваемое на выполнение одного шага расчета в последовательном режиме. Немаловажными параметрами являются также характеристики кластерной системы; в частности, большую роль играет аппаратная структура кластера. В случае, если все используемые в расчете вычислительные модули могут быть размещены на одном узле кластера, итерационная "сшивка" решения требует существенно меньших накладных коммуникационных расходов, что может оказаться критичным для производительности системы, особенно при решении сравнительно небольших задач. На рис. 10 приведен график ускорения вычислении при моделировании обтекания системы круговых цилиндров потоком вязкой несжимаемой среды при значении числа Рейнольдса Не = 1000. Расчетная сетка состояла из 55 000 ячеек.
В пределах одного вычислительного узла (1-4 модуля) ускорение является практически линеиным, а при увеличении числа задеиство-ванных узлов (8, 12, 16 модулеИ) остается почти неизменным. Эффективность использования оборудования существенно повышается при решении более сложных с вычислительнои точки зрения задач. На рис. 11 приведено ускорение вычислении при решении модельнои задачи по расчету конвективного течения воздуха в комнате с подогреваемым полом при проведении расчетов на сетках, содержащих 250 тыс. и 2 млн ячеек. Работа кластера становится близкоИ к оптимальноИ при решении задач, когда на каждыи вычислительныи модуль приходится не менее 30 000 ячеек сетки.
Практикум по параллельным вычислениям. КаждыИ из рассмотренных параллельных вычислительных алгоритмов допускает раз-
4 8 12 16 20 24 28 32 Рис. 11. Графики ускорения расчетов при решении задачи на разных сетках
личные варианты реализации, в том числе с использованием различных технологий параллельного программирования, порой существенно различающиеся по производительности. Более того, эффективность одной и той же реализации во многом определяется архитектурой и конфигурацией конкретного многопроцессорного вычислительного комплекса.
Таким образом, проведение вычислительного эксперимента с использованием высокопроизводительных вычислительных систем требует не только знаний в области архитектуры ЭВМ, умения разрабатывать параллельные алгоритмы, навыков создания параллельных программ, но и способностей планирования и оптимизации загрузки вычислительного комплекса с точки зрения эффективности решения задачи математического моделирования как с позиций минимизации временных затрат, так и в смысле повышения точности и адекватности получаемых результатов. Это означает, что специфика проведения расчетов на многопроцессорных машинах должна приниматься во внимание на всех стадиях математического моделирования [14].
В целях получения студентами базовых знаний в области высокопроизводительных вычислений в учебный план кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана включен Практикум по параллельным вычислениям, первые занятия по которому пройдут в 2014/15 учебном году. В рамках разработки данного курса на кафедре организован семинар "Параллельные вычислительные технологии" для студентов, аспирантов и преподавателей МГТУ им. Н.Э. Баумана, в работе которого также приняли участие представители других вузов и научных организаций (МГУ им. М.В. Ломоносова, Институт системного программирования РАН и др.). Программа семинара включает в себя 11 занятий (табл.2), которые проводятся в компьютерном классе кафедры и сопровождаются практикумом на ЭВМ, в том числе на созданном кластере.
Таблица 2
№ п/п Тема занятия Содержание
1 Введение в параллельные вычисления Актуальные вопросы параллельных вычислений, история и перспективы развития суперкомпьютерной техники, современные технологии разработки масштабируемых приложений, актуальные вычислительно сложные задачи и методы их решения.
2, 3 Особенности языка программирования С++ Обзор языка, базовые и производные типы данных, массивы и указатели, адресация данных в памяти и работа с динамической памятью. Специфика ввода-вывода данных. Функции и аргументы функций. Подходы к разработке переносимых приложений.
4 Архитектура вычислительных комплексов (ВК) Классификация ВК, аппаратная и логическая структура ВК кластерного типа. Классификация Флинна, параллелизм данных и инструкций. Обзор типов вычислителей, многоядерные CPU, видеоускорители GPU. Типы вычислительных сетей, применяемых в ВК. Кластерные операционные системы, оценка производительности кластера.
5, 6 Технология параллельного программирования ОрепМР для систем с общей памятью Обзор технологий параллельного программирования. Многопоточная архитектура параллельной программы. Вычислительный модуль как несколько вычислительных ядер. Основные инструкции библиотеки ОрепМР, спецификаторы переменных, видимость данных и корректность доступа к ним. Методы распараллеливания циклов и контроля распределения загрузки ядер. Примеры параллельных программ: вычисление интеграла, перемножение матриц, решение двумерной нестационарной задачи теплопроводности.
7-9 Технология: параллельного программирования МР1 для систем с распределенной памятью Структура и логика MPI-программы. Основные возможности библиотеки MPI. Типы данных, используемые в MPI, коммуникаторы, группы процессов. Синтаксис основных команд. Обмены данными и синхронизация процессов. Примеры параллельных программ: вычисление интеграла, решение двумерной нестационарной задачи теплопроводности, решение задачи N тел.
10, 11 Технология параллельного программирования CUDA для систем с видеоускорителями nVidia Использование графических ускорителей для вычислений. Особенности архитектуры GPU, специфика разработки алгоритма и организации данных для расчетов на GPU. Библиотека nVidia CUDA, структура CUDA-программы. Синтаксис основных команд, компиляция исходных кодов. Выделение и разработка функции-ядра, исполняемой на GPU, особенности работы с различными типами памяти. Примеры CUDA-программ: различные реализации алгоритма перемножения матриц, решение двумерной нестационарной задачи теплопроводности.
Указанный курс является одним из элементов подготовки специалистов в области математического моделирования, владеющих аппаратом высокопроизводительных вычислений. Дальнейшее развитие
данных способностей у студентов происходит в рамках выполнения курсовых и дипломных работ, нацеленных на решение уже не учебных, а практически значимых задач, связанных с проведением вычислительного эксперимента.
Заключение. Приведены примеры решения практически значимых задач на Учебно-экспериментальном вычислительном кластере кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Проведено исследование эффективности разработанных параллельных алгоритмов и сформулированы рекомендации для повышения производительности кластера. Представлено содержание учебного курса по параллельным вычислениям, ориентированного на студентов кафедры.
Работа выполнена при частичной финансовой поддержке гранта Президента Российской Федерации для государственной поддержки молодых российских ученых — кандидатов наук (проект МК-6482.2012.8), ведущих научных школ Российской Федерации (проекты НШ-1434.2012.2, НШ-255.2012.8), грантов РФФИ (проекты 11-0800699, 12-01-00109, 12-01-31193).
СПИСОК ЛИТЕРАТУРЫ
1. Лукин В. В., Марчевский И. К. Учебно-экспериментальный вычислительный кластер. Часть 1. Инструментарий и возможности // Вестник МГТУ
2. Воеводин В. В., Воеводин В л. В. Параллельные вычисления. СПб.: БХВ-Петербург, 2004. - 608 с.
3. Галанин М. П., Грищенко Е. В., Савенков Е. Б., Токарева С. А. Применение RKDG метода для численного решения задач газовой динамики. Препринт № 52. - М.: ИПМ им. М.В. Келдыша РАН. - 2006. - 31 с.
4. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. -М.: Физматлит, 2001. - 608 с.
5. MignoneA., Bodo G. Shock-capturing schemes in computational MHD //
6. TorrilhonM. Locally divergence-preserving upwind finite volume schemes for
7. Марчевский И. К., Токарева С. А. Сравнение эффективности параллельных алгоритмов решения задач газовой динамики на разных вычислительных комплексах // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки.
8. Т о t h G. The div В = 0 constraint in shock-capturing magnetohydrodynamics
9. Галанин М. П., Лукин В. В. Разностная схема для решения двумерных задач идеальной МГД на неструктурированных сетках. Препринт № 50. - М.: ИПМ им. М.В. Келдыша РАН. - 2007. - 29 с.
10. CottetG. -H., Koumoutsakos P. D. Vortex methods: Theory and practice.
11. ДынниковаГ. Я. Вихревые методы исследования нестационарных течений
вязкой несжимаемой жидкости: Дис. ... д-ра физ.-мат. наук. - М., 2011. - 269 с.
12. М а р ч е в с к и И И. К., М о р е в а В. С. ПараллельныИ программный комплекс РОЬАКА для моделирования обтекания профилеИ и исследования расчетных схем метода вихревых элементов // Параллельные вычислительные технологии (ПаВТ-2012): Труды международноИ научноИ конференции (Новоси-
13. МарчевскиИ И. К., Щеглов Г. А. Модель симметричного вортона-отрезка для численного моделирования пространственных течениИ идеальноИ несжимаемоИ среды // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные на-
14. 3 а р у б и н В. С. Математическое моделирование в технике. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2001. - 496 с.
Статья поступила в редакцию 25.06.2012
