CC BY
63
ISSN 1992-6502 ( P ri nt)_
2014. Т. 18, № 2 (63). С. 218-226
Ъыьмт QjrAQnQj
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 519.6
Параллельный алгоритм RKDG метода второго порядка
для решения двумерных уравнений идеальной магнитной гидродинамики
1 2 м. п. Галанин , в. в. Лукин , к. л. Шаповалов
1 [email protected], 2 [email protected]
1,2 ФГБУН «Институт прикладной математики им. М. В. Келдыша РАН» 1,2 ФГБОУ ВПО «Московский государственный технический университет им. Н. Э. Баумана»
Поступила в редакцию 22.12.2013
Аннотация. Рассмотрена параллельная реализация RKDG метода решения двумерных уравнений идеальной магнитной гидродинамики на треугольных неструктурированных сетках. Описан программный комплекс, построенный с применением технологий MPI и OpenMP. Приведены алгоритмы монотонизации решения и бездивергентной реконструкции магнитного поля, позволяющие получать физически адекватные результаты расчетов с высоким порядком точности. Обсуждены результаты тестовых расчетов и полученные значения ускорения вычислений за счет использования параллельных технологий.
Ключевые слова: магнитная гидродинамика; RKDG метод; параллельное программирование; MPI; OpenMP.
Уравнения идеальной магнитной гидродинамики (МГД) описывают множество явлений, представляющих как чисто научный, так и прикладной технический интерес. В частности, в рамках моделей, основанных на системе уравнений МГД, рассматриваются процессы электромагнитного ускорения вещества в различных условиях (плазменные ускорители, астрофизические джеты). Особый интерес представляет моделирование развития неустойчивости течения замагниченной плазмы в окрестности поверхности разрыва (ударной волны), для адекватного описания которого целесообразно использовать численные методы второго и более высоких порядков точности.
Одним из наиболее современных и популярных при решении задач газовой динамики и МГД является метод RKDG [1, 2] (Runge-Kutta Discontinuous Galerkin). Этот метод обеспечивает высокий уровень разрешения разрывов и позволяет повышать порядок точности метода без увеличения шаблона аппроксимации, что осо-
Выполнено при финансовой поддержке РФФИ (проекты 12-01-00109, 14-01-31496, 12-02-12096, 12-0200687), гранта Президента РФ для государственной поддержки ведущих научных школ РФ (проект НШ-6061.2014.2).
Рекомендовано к публикации Программным комитетом международной конференции ПаВТ'2013.
бенно важно при использовании неструктурированных сеток. В то же время методы высокого порядка являются немонотонными и могут приводить к накоплению численного магнитного заряда, что создает существенные сложности при моделировании замагниченных течений.
В работе рассмотрена модификация разрывного метода Галеркина второго порядка для решения системы уравнений двумерной магнитной гидродинамики на треугольных сетках. Построенный метод включает процедуры монотонизации численного решения, а также позволяет получать решение, удовлетворяющее условию бездивергентности магнитного поля.
Анализ процесса формирования неустойчивости течения замагниченной плазмы требует использования подробных сеток. В то же время применение RKDG метода само по себе является вычислительно затратным и приводит к необходимости использования параллельных вычислительных технологий для систем с общей и распределенной памятью [3, 4]. В работе описана параллельная реализация данного метода. Показана эффективность распараллеливания на тестовых задачах МГД.
СИСТЕМА НЕСТАЦИОНАРНЫХ
УРАВНЕНИЙ ИДЕАЛЬНОЙ МАГНИТНОЙ ГИДРОДИНАМИКИ
Система нестационарных уравнений идеальной магнитной гидродинамики, записанная в консервативных переменных, включает в себя законы сохранения массы, импульса, энергии, а также уравнение Фарадея для магнитной индукции [5]:
где
ди дРл дР?
и = (р, руг, ру2, ру з, е, Въ В2, В3)т,
(1)
В'
В(
= (рк 1, ру1 + Р + -- тЬ'Р^Ъ -
в1в3
В, в
1°2
8п 4п В2
рргр з
4ттг ' (В-у)
-—, (е + р + —)гЛ| — Вл —--,
477: р 8тГ 1 1 4л:
О, ргВ2 ~ у2В1, угВ3 - У^У,
В?В-1
В'
В4
Р2 = (р*2, р^ - _ р*2 + р + - - - ,
В2В3 В2 (В • V)
у2В1 - г^Вг, 0, У2В3 - р3В2)т.
Здесь р - плотность газа, V = (V}, у2 , у3 ) т -вектор скорости, р - давление, е - полная энергия единицы объема газа, В = (Вг, В2, В3 )т -вектор магнитной индукции. Из закона Фарадея следует условие бездивергентности магнитного поля с1 IV В = 0. Будем считать, что плазма является сильно ионизированным, неполяри-зованным совершенным газом с уравнением состояния , где - показатель
адиабаты, £ - удельная внутренняя энергия. При
Р , ру2 . в2
этом е =--\---\--.
у-1 2 8п
Система (1) с уравнением состояния газа является замкнутой. Можно показать, что эта система является гиперболической [5] и, следовательно, существует разложение
ОД
—Г = А; = ди 1 1 1 1
I = 1,2,
(2)
где и - ортогональные матрицы,
составленные из левых и правых собственных векторов, а - действительная диагональная матрица, составленная из собственных чисел матрицы .
ЧИСЛЕННЫЙ МЕТОД
Разрывный метод Галеркина. Для
решения системы (1), или в индексной форме записи (здесь и далее в этом параграфе принято
правило суммирования по повторяющимся индексам)
Зи,- ОД-; — + = 0, дЬ Зху
(3)
на треугольной сетке используем разрывный метод Галеркина [1, 2]. Применим стандартную схему методов галеркинского типа - умножим правую и левую часть (1) на произвольную скалярную функцию , проинтегрируем по объему , применим правило интегрирования по частям и получим:
_ Г р ^М = 0, (1) (4)
■^д ]>1 дХ] ' к '
где - единичный вектор внешней нормали к границе области - поверхности .
Рассмотрим двумерную область И и соответствующую ей трехмерную область Ъ = й X М. Введем в области й сетку П состоящую из системы треугольников
. В качестве области рассмотрим единичную треугольную призму с сечением . Тогда второй интеграл в (4) в силу свойств системы (1) можно преобразовать к виду [6]
/
,, иа
л
= X / ^"" ТкМи1' ■■ )ШС1Б'
к=1 сг
где - поверхность -й боковой грани призмы , - матрица, осуществляющая
поворот векторов и на угол, равный углу между нормалью и положительным
направлением оси хъ а Р* - вектор, описывающий поток консервативных переменных через поверхность . Будем искать решение в виде ^ лгй
и = ^ ^ Уа,р(р')<ра1р(х1,х2'),
а=1 у6=1
где - зависящее от порядка аппроксимации число базисных функций, определенных в ячейке К а; { р а,р (х1 ,х{) х -ортонормированная система базисных функций (полиномов степени не выше ),
определенных в а-й ячейке; уа'@ (Ь) - вектор коэффициентов при базисных функциях. В качестве пробных функций выберем те же
базисные функции ра. р. Тогда приходим к системе обыкновенных дифференциальных уравнений разрывного метода Галеркина порядка
г
+ / г^-о. (5)
уа '
а = Т^Щ, р = I = 1Д
где хр - номер ребра а-й ячейки, / (а. хр)- номер соседней с ней ячейки, Г* (уа. у1 ) -численный поток, зависящий от значений решения в ячейках, примыкающих к ребру , и определяемый из решения соответствующей задачи Римана о распаде разрыва [5].
Система (5) дополняется начальным условием вида
Уа,р(о ) - | ЩРа,р<М. (6)
Vя
где - начальное распределение
консервативных переменных.
Интегралы
по поверхностям и по
областям в рассматриваемом двумерном случае сводятся к интегралам по ребрам и треугольнику Ка соответственно.
Для приближенного решения задачи Римана для системы уравнений МГД используется ряд численных методов, таких как HLL, HLLC, ИЬЬБ [5, 7]. Метод HLLD позволяет получить наилучшее разрешение разрывов, но во многих случаях приводит к возникновению численной немонотонности в решении. Чтобы этого избежать, может быть использован численный поток HLL, обладающий большей численной вязкостью.
При расчете интегралы заменяются квадратурными формулами, точность которых согласована с порядком метода . Решение задачи Коши (5)-(6) производится численным интегрированием явным методом Рунге-Кутты, шаг по времени определяется динамически из условия устойчивости Куранта-Фридрихса-Леви:
т ах | Л( <1 . (7)
1 "тт
где - длина наименьшего ребра ячейки в сетке, - собственные числа системы (1) на предыдущем шаге по времени. Разрывный метод Галеркина первого порядка совпадает с соответствующими методами типа Годунова [5].
Для аппроксимации граничных условий используется метод фиктивной ячейки.
Монотонизация решения. Решение уравнений МГД разрывным методом Галеркина второго порядка точности требует использования функций-ограничителей наклона решения в ячейке - лимитеров [6] - как для поддержания монотонности решения, так и для предотвращения появления на очередном временном слое нефизичных отрицательных значений плотности и давления.
В данной работе использован «ТУВ minmod» лимитер для кусочно-линейных функций [2, 5], значение которого зависит от исходного наклона решения в данной ячейке и от средних значений решения в соседних ячейках. Процедура лимитирования магнитного поля может привести к возникновению численного магнитного заряда, и поэтому применяется только при вычислении потоков. Лимитирование проводится для локальных характеристических переменных, получаемых из консервативных умножением на матрицу левых собственных векторов, вычисляемую в каждой ячейке вместе с разложением (2).
Рассмотрим алгоритм применения лимитера к вектору консервативных переменных, заданному линейной функцией на
треугольнике К 0 с соседями К 1 , К2, К3 (рис. 1).
Рис. 1. Построение лимитера
Для каждой пары соседних треугольников , (рис. 2) из системы уравнений П1( - Ь0 = ачД(Ь( - Ь0) + «¿у,2(Ь/ - Ь0), где - барицентр -го треугольника, -центр -го ребра, общего между и , определим скалярные коэффициенты и
, характеризующие геометрические свойства данной группы треугольников , , .
На основании полученных коэффицентов определим средние значения скачка функции ико в т (:
и
о(Ш() = ai}д (%i(b() - %o(b0)) +
+aijt2(uKj(<by) - uKo(b0)).
Рис. 2. Схема перераспределения магнитных потоков
Исходя из полученных значений, проведем ограничение скачков в центрах ребер:
Д( = Кт(1(ико(Ш() -
—ико (Ь0)), уЬЙ^о (ГП())< * = 1,3,
где - константа, определяющая степень
сглаживания решения, и - матрицы левых и правых собственных векторов из разложения (2), т - «ТУВ minmod» функция, применяемая покомпонентно для каждой характеристической переменной.
т
lj(ä,b) =
iy,если|ау| < Mj(hminy,
0.5 i sign clj + sign bj ) min(| aj |, |by |),
7 = 1,8,
где - константа, зависящая от задачи, - функция знака. После определения скачков , задающих новый наклон решения в , их значения используем для обновления коэффициентов при базисных функциях.
2.3. Условие отсутствия магнитного заряда
Выполнение условия бездивергентности магнитного поля является одной из
основных проблем при создании численного метода решения уравнений магнитной гидродинамики [8]. Несмотря на то, что дивергенция начального распределения магнитного поля равна нулю, при использовании консервативной формы записи системы уравнений МГД (1) в ходе расчетов может накапливаться численный магнитный заряд. Такая ситуация приводит к нефизичным результатам или даже прежде-
временному прерыванию расчета из-за получения отрицательных значений плотности или давления.
Для поддержания бездивергентности магнитного поля нами реализован метод, основанный на использовании -компоненты электрического поля и ее градиента в узлах сетки (рис. 2) для определения бездивергентного потока магнитного поля. Для случая метода первого порядка (без использования градиентов) процедура описана в [9].
После определения известными методами (HLL, HLLC, HLLD) значений численного потока в точках интегрирования (для
метода второго порядка две точки - )
-компоненты электрического поля и
_ /дБ? дЕЛ. У с 2 - V аж/ажгЛ
ik
+
k=l
+
а, \ (8)
(F6*[ife,2],F7*[^2])-eifc
I
k=1
Здесь и далее е ^ - единичный вектор, касательный к ребру и направленный из -й вершины в -ю, - длина ребра , - число ребер, выходящих из вершины . Восстановление градиента производится по фор-1 1 муле « 1 /5 - 1 Я2п С / (рис. 2).
Этот интеграл берется численно с использованием вычисленных значений потоков
Р* [ ц,к ] .
Для монотонизации решения необходимо провести ограничение градиента электрического поля аналогично процедуре, рассмотренной выше. Для этого на основании значений в узлах расчитываем средний градиент , результирующее значение
градиента, используемое далее в расчетах, определяем путем применения minmod функции .
Используя значения электрического поля и его градиента, можно построить бездивергентные потоки и , отличающиеся от исходного только компонентами потока магнитного поля:
ТЦЦ, к] = к], к = 1,2,1 = 1,5,8, {П[ЧЛ1Т*7[Ц,1])Т =
= - ((3 + л/3)Ея + (3 - л/3)Ег] + 1
+ 2 - ' е(у/г(у) кце1},
= g((3-V3)£zi + (3 + V3)£zy + 1
+ 2 (VEzi - VEzy) • etjhtj) htjetj,
Такое перераспределение потоков сохраняет консервативность схемы и гарантирует без-дивергентность решения. Это подтверждается численными экспериментами: в случае метода первого порядка этот факт показан в [9]. В то же время численные эксперименты показывают, что такой подход приводит к сглаживанию разрывов (вносит в схему дополнительную численную взякость): на 4-5 ячеек для разрывного метода Галеркина первого порядка, и 2-3 ячейки для метода второго порядка.
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ
Описанный численный метод реализован в параллельном программном комплексе для кластерных систем с разделенной памятью. Комплекс написан на языке Fortran с использованием технологий OpenMP и MPI. Параллельный алгоритм комплекса основан на идее разбиения исходной расчетной области на подобласти, решение в каждой из которых рассчитывается отдельным вычислительным модулем [10]. Поскольку численный метод -явный, подобный подход позволяет добиться высокой эффективности распараллеливания. Целостность решения обеспечена наличием обменов смежными граничными данными между модулями, обмены реализованы средствами технологии MPI. Расчет внутри каждой из подобластей может быть ускорен с помощью применения технологии OpenMP в том случае, если задействованные для решения задачи модули состоят из нескольких вычислительных ядер [3].
Опишем основные этапы работы параллельного программного комплекса.
Подготовка к запуску. Расчет производится на треугольных неструктурированных сетках. Для построения сеток использован программный комплекс Gridder2D [11], позволяющий производить триангуляцию областей сложной формы. На подготовительном этапе производится разбиение сетки на подобласти по числу используемых вычислительных модулей. При этом необходимо, с одной стороны, минимизировать количество смежных данных, подлежащих обмену между модулями, а с другой - сбалансировать число ячеек сетки в каждой подобласти для обеспечения равномерной загрузки ядер. Эту задачу можно представить как построение
разбиения графа, каждая вершина которого соответствует ячейке, а ребро между вершинами - общему ребру между ячейками. Для приближенного решения данной задачи используется программа Metis [12], устанавливающая соответствие между номером ячейки в сетке и номером подобласти.
Инициализация данных. Инициализация расчета производится процессом, исполняемым на нулевом модуле («корневой процесс»). На основании информации о расчетной сетке и ее разбиении на подобласти, а также о сеточном шаблоне, используемом как для расчета решения на очередном временном слое, так и для его лимитирования, корневой процесс формирует для каждого расчетного модуля ряд вспомогательных массивов.
• Массив перенумерации узлов, устанавливающий соответствие между номерами узлов в исходной сетке и номерами узлов, используемых в модуле - на каждом вычислительном модуле хранится и участвует в расчетах лишь минимально необходимая часть сетки и определенные на ней распределения физических переменных задачи.
• Массив перенумерации ячеек сетки -аналогично массиву перенумерации узлов данный массив позволяет хранить в памяти информацию только о тех ячейках, которые необходимы для расчета на данном модуле. Этот массив специальным образом упорядочен. Самые маленькие номера имеют ячейки, которые используются только на данном модуле. Далее размещены блоки ячеек, значение решения в которых рассчитывается данным модулем и отсылается на другой модуль путем групповой пересылки. Затем располагаются ячейки, значение решения в которых по отдельности рассылается нескольким модулям (угловые ячейки). В конце массива расположены блоки ячеек, значения решения в которых принимаются от других модулей блочным способом, ячейки, решение в которых принимается по отдельности, а также граничные ячейки. Такое упорядочивание позволяет достичь большей эффективности пересылки данных между модулями.
• Массивы, описывающие начало и конец рассылаемых и принимаемых блоков и одиночных данных о решении.
После рассылки данных корневым модулем все вычислительные модули (включая корневой) начинают независимый расчет решения в своей подобласти, инициализируя массивы, описывающие базисные функции, содержащие
скалярные коэффициенты лимитера и ряд других.
Для обмена значениями решения в смежных ячейках и электрического поля в смежных узлах на границах подобластей каждым модулем инициализируются отложенные передача и прием данных. При этом благодаря непрерывной блочной структуре хранения данных в памяти пересылка значений в смежных ячейках и узлах осуществляется без формирования дополнительных буферных массивов (и без дополнительных накладных расходов) прямым обращением к соответствующим указателям.
Так как z-компонента электрического поля, ее производные и значения лимитера для этих производных вычисляются в узлах сетки по формуле (8) как сумма значений на ребрах, прилегающих к данному узлу, то при параллельном счете требуется учитывать значения этих величин на ребрах, обрабатываемых разными модулями. В этих случаях каждый из модулей вычисляет сумму только для своей части ребер и рассылает отложенным способом это значение на другие модули, обрабатывающие данный узел. Затем каждый из модулей суммирует принятые данные и получает правильные значения в узле. Такая методика оказывается эффективнее централи-зованого расчета соответствующих значений.
Расчет. Интегрирование по времени производится явным методом Рунге-Кутты второго порядка: каждый шаг по времени состоит из двух этапов вычисления решения в моменты времени, соответствующие выбранному методу Рунге-Кутты. Величина очередного шага по времени для метода Рунге-Кутты в каждой подобласти определяется из условия (7), а из этих шагов при помощи средств библиотеки MPI выбирается минимальный, при котором условие устойчивости численной схемы оказывается выполнено во всей расчетной области. Процедура расчета решения в очередной временной точке включает в себя как пересчет решения внутри ячеек подобласти, так и вычисление лимитеров для всех консервативных переменных, а также операции по вычислению в узлах сетки электрического поля, его градиента и соответствующих ограничителей наклона. Все указанные действия требуют обмена данными между смежными вычислительными модулями, которые выполняются в режиме отложенных запросов, что позволяет избежать простоя вычислителей. Запись результатов расчета на диск производится буфери-
зованным способом. В случае, когда вычислительный модуль включает несколько вычислительных ядер, последние два приема позволяют существенно уменьшить реальную долю последовательного кода во времени исполнения алгоритма. Кроме того, расчеты, выполняемые внутри одного вычислительного модуля, в этом случае легко распараллелены при помощи технологии OpenMP.
РЕЗУЛЬТАТЫ РАСЧЕТОВ
Созданный программный комплекс протестирован на ряде стандартных тестовых задач двумерной идеальной МГД. Результаты решения задачи о вихре Орзага-Танга, демонстрирующие разрешающую способность метода и эффективность распараллеливания алгоритма, приведены ниже. Результаты получены RKDG методом второго порядка точности по пространству и времени с HLLD разностными потоками.
Вихрь Орзага-Танга. Данная задача иллюстрирует процесс формирования и взаимодействия ударных волн [13]. Начальные распределения физических величин - периодические по пространству гладкие функции. Тем не менее достаточно быстро в решении появляются разрывы плотности, давления и скорости. Благодаря периодичности решения по пространству в качестве расчетной области достаточно взять единичный квадрат с периодическими граничными условиями. Использованы следующие начальные данные
Go,v1,v2,v3,p,B1,B2,B3) = 25
= —sin277:x2, sin277:x1, О,
36л:
5
-, — sm2nx7, sin4nx,, 0).
12л: 2 1 ^
Показатель адиабаты , что соот-
ветствует идеальному одноатомному газу. Число Куранта выбрано равным С = 0 . 4, в расчетах использовались треугольные сетки различной подробности. На рис. 3 представлены распределения плотности и модуля вектора магнитной индукции в момент времени t = 0 . 5 , полученные на сетке из 903222 треугольников. Метод показывает высокий уровень разрешения структуры решения, не сравнимый с методами первого порядка (см. аналогичные результаты в [13]), причем это касается не только обострения разрывов (в среднем разрывы размазаны на 2-3 ячейки), но и выявления мелкомасштабной структуры разрывов в решении.
Рис. 3. Распределение плотности (слева) и модуля вектора магнитной индукции (справа)
в задаче о вихре Орзага-Танга
Параллельные расчеты. Проведено тестирование параллельного алгоритма на гибридном вычислительном кластере К-100 Института прикладной математики им. М.В. Келдыша РАН. Оценка ускорения проводилась для решения на сетках из 585508 и 903222 элементов. Доля параллельного кода в описанном алгоритме при решении задачи о вихре Орзага-Танга составляет 99.98 % для сетки из 585508 элементов и 99.99 % - для 903222 элементов. Согласно закону Амдала [14], такие доли ограничивают максимальное ускорение расчетов на 256 ядрах относительно последовательной версии величинами 243 и 249 раз соответственно. Решение задачи на сетке из 903222 элементов на одном ядре заняло 1843,8 мин, на 256 ядрах - 8,17 мин (ускорение в 225.7 раза). На сетке из 585508 элементов на одном ядре - 720 мин, на 256 ядрах - 3,68 мин (ускорение в 195.7 раза).
Расхождение реального ускорения с оценкой по закону Амдала объясняется главным
образом не влиянием межмодульных обменов, а спецификой аппаратной структуры кластера. Вычислительные модули кластера К-100 сгруппированы по 11 на вычислительный узел (в одноядерной конфигурации модуля), каждый из которых представляет собой сервер с двумя шестиядерными процессорами Intel Xeon 5670. Подобная плотность вычислительных модулей в расчете на узел приводит к конкуренции запросов модулей одного узла к оперативной памяти, а следовательно - к большому числу кэш-промахов. Это особенно критично для текущей реализации метода решения уравнений МГД второго порядка точности, поскольку каждой ячейке сетки соответствует как минимум 32 числа с плавающей точкой. Расчеты показали, что процедура вычисления потоков через ребра на четырех ядрах ускоряется в 3,8 раза, в то время как процедура копирования решения с предыдущего слоя на новый - только в 2,1 раза.
чшш
■-г-•-- ■ , 1 ■ , т -,-,-,-
1 2 4 3 16 32 64 128 255
* Ускорение расчетов на сетке 903 тысячи элементов 1 ■ Ускорение расчетов на сетке 535 тысяч элементов Линейное ускорение • Ускорение расчетов по закону Амдалу
Рис. 4. Ускорение расчетов в задаче о вихре Орзага-Танга на различных сетках
0,7 -1-1-1-1-1-1-1-1-
1 2 4 8 16 32 64 12S 2S6
• Эффективность на сетке 903 тысячи элементов ■ Эффективность на сетке 585 тысяч элементов —А—Эффективность прилинейном ускорении • Эффективность по закону Амдэлу
Рис. 5. Эффективность параллельного алгоритма в задаче о вихре Орзага-Танга на различных сетках
Специфика аппаратной структуры кластера К-100 объясняет и наблюдаемое на рис. 4 и 5 сверхлинейное ускорение расчетов, возникающие на 2 и 4 ядрах. В этом случае не возникает конкуренции за ресурсы между вычислительными модулями одного узла, причем часть обменов может выполняться незанятыми в расчете ресурсами узла.
На рис. 4 представлены графики ускорения расчетов, полученные при решении задачи о вихре Орзага-Танга на сетках из 585508 и 903222 элементов. Также представлен график ускорения согласно закону Амдала для 99.98 % параллельного кода и график линейного ускорения. На рис. 5 представлены графики эффективности (отношение ускорения к числу задействованных ядер) распараллеливания расчетов на различных сетках, эффективности, согласно закону Амдала, для 99.98 % параллельного кода и эффективности при линейном ускорении. Видно, что эффективность распараллеливания при большом числе используемых вычислительных модулей существенно зависит от трудоемкости задачи. В целом метод показывает высокую эффективность распараллеливания и хорошую масштабируемость вплоть до сотен задействованных вычислительных модулей.
ЗАКЛЮЧЕНИЕ
Создан алгоритм для решения двумерных уравнений магнитной гидродинамики разрывным методом Галеркина, основывающийся на методах приближенного решения задачи Римана о распаде разрыва (HLL, HLLC, HLLD). Рассмотрены способы монотонизации схемы, введена процедура реконструкции магнитных потоков, позволяющая в случае методов первого и второго порядков точности получать бездивергентные распределения магнитного
поля. Протестирован метод второго порядка точности по пространству и времени. Тестовая задача о вихре Орзага-Танга показывает высокий уровень разрешения разрывов и монотонность получаемого численного решения. Создана и протестирована параллельная версия программы с использованием технологии параллельных вычислений для систем с распределенной памятью MPI. Представлены графики ускорения и эффективности распараллеливания расчетов для тестовой задачи, показывающие высокую эффективность и хорошую масштабируемость параллельного алгоритма.
СПИСОК ЛИТЕРАТУРЫ
1. Cockburn B., Shu C.-W. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems // Journal of Scientific Computing. 2001. № 3 (16). P. 173-261. [ B. Cockburn and C. W. Shu, "Runge-Kutta discontinuous Galerkin methods for convection-dominated problems," Journal of Scientific Computing, no. 3 (16), pp. 173-261, 2011. ]
2. Галанин М. П., Савенков Е. Б., Токарева С. А. Решение задач газовой динамики с ударными волнами RKDG-методом // Математическое моделирование. 2008. Т. 20, № 11. C. 55-66. [ M. P. Galanin, E. B. Savenkov, and S. A. Tokareva, "The solution of gas dynamics problems with shock waves using Runge-Kutta discontinous Galerkin method," (in Russian), Matematicheskoe modelirovanie, vol. 20, no. 11, pp. 55-66, 2008. ]
3. Лукин В. В., Марчевский И. К. Учебно-экспериментальный вычислительный кластер. Часть 1. Инструментарий и возможности // Вестник МГТУ им. Н. Э. Баумана. Сер. Естественные науки. 2011. № 4. С. 28-44. [ V. V. Lukin and I. K. Marchevskii, "Computing Cluster for Training and Experiments. Part 1. Tools and Possibilities", (in Russian), Vestnik MGTU im. N. E. Baumana, Estestvennye nauki, no. 4, pp. 28-44, 2011. ]
4. Учебно-экспериментальный вычислительный кластер. Часть 2. Примеры решения задач / В. В. Лукин [и др.] // Вестник МГТУ им. Н. Э. Баумана. Естественные науки. 2012. № 4. С. 82-102. [ V. V. Lukin, et al., "Computing Cluster for Training and Experiments. Part 2. Examples of Solving
Problems", (in Russian), Vestnik MGTU im. N. E. Baumana, Estestvennye nauki, no. 4, pp. 82-102, 2012. ]
5. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю.
Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с. [ A. G. Kulikovskii, N. V. Pogorelov, and A. Yu. Semenov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, (in Russian). Moscow: Fizmatlit, 2001. ]
6. Toro E. F. Riemann solvers and numerical methods for fluid dynamics. A practical introduction. Berlin: Springer, 2009. 724 pp. [ E. F. Toro, Riemann solvers and numerical methods for fluid dynamics. A practical introduction. Berlin: Springer, 2009. ]
7. Mignone A., Bodo G. Shock-capturing schemes in computational MHD // Lect. Notes Phys. 2008. V. 754. P. 71101. [ A. Mignone and G. Bodo, "Shock-Capturing Schemes in Computational MHD," Lect. Notes Phys., vol. 754, pp. 71101, 2008. ]
8. Toth G. The div B = 0 constraint in shock-capturing magnetohydrodynamics codes // J. of Comp. Physics. 161, 2000, p. 605-652. [ G. Toth, "The div B = 0 constraint in shock-capturing magnetohydrodynamics codes," J. of Comp. Physics, no. 161, pp. 605-652, 2000. ]
9. Torrilhon M. Locally divergence-preserving upwind finite volume schemes for magnetohydrodynamic equations // J. Sci. Comp. 2005. Vol. 26. P. 1166-1191. [ M. Torrilhon, "Locally Divergence-preserving Upwind Finite Volume Schemes for Magnetohydrodynamic Equations," J. Sci. Comp., vol. 26, pp. 1166-1191, 2005. ]
10. Марчевский И. К., Токарева С. А. Сравнение эффективности параллельных алгоритмов решения задач газовой динамики на разных вычислительных комплексах // Вестник МГТУ им. Н. Э. Баумана. Естественные науки. 2009. № 1. С. 90-97. [ I. K. Marchevskii and S. A. Tokareva, "Comparison of efficiency of parallel algorithms to solve problems of gas dynamics at various computing complexes," (in Russian), Vestnik MGTU im. N.E. Baumana, Estestvennye nauki, no. 1, pp. 90-97, 2009. ]
11. Щеглов И. А. Программа для триангуляции сложных двумерных областей Gridder2D: препр. Ин-та прикл. матем. им. М. В. Келдыша РАН. 2008. № 60. 32 с. [ I. A. Sheglov, Software for constrained 2D triangulation Gridder2D, Keldysh Institute preprints, no. 60, 2008. ]
12. Karypis G., Kumar V. A fast and high quality multilevel scheme for partitioning irregular graphs // SIAM Journal on Scientific Computing. 1999. V. 20, no. 1. P. 359-392. [ G. Karypis and V. Kumar, "A fast and high quality multilevel scheme for partitioning irregular graphs," SIAM Journal on Scientific Computing, vol.. 20, no. 1, pp. 359-392, 1999. ]
13. Галанин М. П., Лукин. В. В. Разностная схема для решения двумерных задач идеальной МГД на неструктурированных сетках: препр. Ин-т прикл. матем. им. М. В. Келдыша РАН, 2007. № 50. 29 с. [ M. P. Galanin and V. V. Lukin, Finite-difference scheme for the two-dimensional ideal MHD problems using unstructured grids, Keldysh Institute preprints, 2007, no. 50. ]
14. Воеводин В. В., Воеводин Вл. В. Параллельные вычисления. СПб.: БХВ-Петербург, 2004. 608 с. [ V. V. Voevodin and Vl. V. Voevodin, Parallel computing, (in Russian). Saint-Petersburg: BHV-Peterbuirg, 2004. ]
ОБ АВТОРАХ
ГАЛАНИН Михаил Павлович, зав. отделом, проф. каф. прикладной математики. Дипл. физик (МГУ им. М. В. Ломоносова, 1979). Д-р физ.-мат. наук по матем. моделированию (там же, 1995). Иссл. в обл. матем. моделирования и числ. методов.
ЛУКИН Владимир Владимирович, науч. сотр. Дипл. инж.-математик (МГТУ им. Н. Э. Баумана, 2008). Канд. физ.-мат. наук по матем. моделированию (ИПМ им. М. В. Келдыша, 2011). Иссл. в обл. матем. моделирования и числ. методов.
ШАПОВАЛОВ Кирилл Ладиславович, дипл. инж.-математик (МГТУ им. Н. Э. Баумана, 2013). Иссл. в обл. числ. методов магнит. гидродинамики.
METADATA
Title: Parallel algorithm for second order RKDG method for 2D
magnetohydrodynamics system of equations. Authors: M. P. Galanin, V. V. Lukin, and K. L. Shapovalov. Affiliation:
Keldysh Institute of Applied Mathematics RAS, Russia. Bauman Moscow State Technical University (MGTU), Russia. Email: [email protected]. Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 18, no. 2 (63), pp. 218-226, 2014. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: A parallel implementation of RKDG method for solving two-dimensional ideal MHD equations on triangular unstructured grid is considered. A software package using technologies MPI and OpenMP is described. The algorithms for solution monotonization and divergence-free reconstruction of the magnetic field for physically adequate results of calculations with a high order accuracy obtaining are presented. The results of test calculations acceleration obtained through the use of parallel technologies are discussed. Key words: magnetohydrodynamics; RKDG method; parallel programming; MPI; OpenMP.
About authors:
GALANIN, Michail Pavlovich, Head of Dept., KIAM RAS, Prof., Dept. of Applied Math. Dipl. Physicist (Lomonosov MSU, 1979). Dr. of Phys.-Math. Sci. (Lomonosov MSU, 1995).
LUKIN, Vladimir Vladimirovich, Researcher. KIAM RAS, Assoc. Prof, Dept. of Applied Math. Dipl. Engineer-math. (MGTU, 2008). Cand. of Phys.-Math. Sci. (KIAM RAS, 2011).
SHAPOVALOV, Kirill Ladislavovich, Dipl. Engineer-math. (MGTU, 2013).
