CC BY
70
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
УДК 524.6+524.7
С. А. Хоперсков, А. В. Хоперсков, А. В. Засов, М. А. Бутенко
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ МОДЕЛИРОВАНИЯ ДИНАМИКИ ГАЗА В СИЛЬНО НЕОДНОРОДНЫХ ГРАВИТАЦИОННЫХ ПОЛЯХ
Создан пакет программ, реализующий трехмерную численную схему TVD типа MUSCL для решения уравнений газовой динамики во внешних гравитационных полях с учетом тепловых процессов и самогравитации в декартовой и цилиндрической координатных системах. Численный алгоритм адаптирован для расчетов на компьютерах с массивно-параллельной архитектурой. Основываясь на построенной численной модели, исследованы физические процессы, протекающие в галактических газовых дисках. Детально изучен механизм гофрировочной неустойчивости спиральной структуры в газовом диске, рассмотрен новый гидродинамический механизм образования кольцеобразных галактических структур. Проведено моделирование облачной структуры газопылевой подсистемы Млечного Пути, показавшее хорошее согласие с данными наблюдений. Газодинамика; численные методы; параллельные алгоритмы; физика галактик
Изучение различных аспектов феномена спиральных галактических структур сохраняет свою актуальность, несмотря на 50-летнюю историю исследований [1]. В основе построения современных моделей спиральных галактик лежат различные методики компьютерного моделирования. Заметный прогресс был достигнут с привлечением высокопроизводительных суперкомпьютеров и использованием параллельных технологий.
Физические параметры, характеризующие свойства галактического газа, лежат в очень широких пределах, и не имеют аналогов при решении типичных задач, возникающих в технике или «земной» физике.
Межзвездная среда, являясь многофазной, имеет температуру от 3 К до 105 К и перепады плотности, различающиеся на несколько порядков [2]. Типичные значения скорости составляют ~ 100-300 км/с, что дает числа Маха 10 ~ ~ 100. Наблюдается иерархия пространственных и временных масштабов, которые существенно определяют свойства галактической системы. Например, газовый диск крупной спиральной галактики простирается до 20-40 кпк по радиусу. С другой стороны, мелкомасштабные структуры, формирующиеся в области спиральных
Контактная информация: [email protected] Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2011». Работа была выполнена при финансовой поддержке грантов РФФИ № 10-07-97017, гранта ВолГУ № 70-2011-а/ВолГУ, некоммерческого фонда «Династия», ФЦП «Старт-11 Н1» (проект № 8861р/14383) и ФЦП Росбразование госконтракт П1248 (ФЦП НК).
галактических ударных волн, составляют всего ~ 10 пк, и для их качественного разрешения необходимы масштабы до ~ 1 пк.
По современным представлениям галактики являются существенно нестационарными объектами [3]. Характерные времена процессов, определяющих динамику газового диска в целом, достигают 1 млрд лет, а необходимость одновременного моделирования быстропроте-кающих тепловых явлений, нелинейных этапов развития физических неустойчивостей (гравитационных и сдвиговых) требует рассмотрения динамики на временах до 104 лет.
Отличительной особенностью моделирования динамики межзвездной среды в пределах галактики является богатый спектр физических эффектов и процессов, существенно влияющих на динамику газа. Отметим целый ряд газодинамических и плазменных неустойчивостей различной природы: гравитационная, тепловая, конвективная, магнитная. Это вызывает необходимость дополнительного учета в уравнениях гидродинамики тепловых процессов, химических превращений, самогравитации, что делает задачу очень жесткой. Отмеченные выше временные характеристики задачи требуют большого числа временных шагов интегрирования (> 105). Например, для самосогласованного описания динамики звездно-газового диска на современном научном уровне, необходимо для моделирования звездной подсистемы более 107 гравитационно взаимодействующих частиц [4], а число ячеек сетки для газовой компоненты должно превышать 108, что требует значительных ресурсов памяти и машинного времени. Без использования параллельных технологий реше-
ние такого рода задач оказывается практически невозможным.
Выделим основные направления работы по исследованию динамики галактик с использованием параллельных технологий:
• создание пакета программ для численного моделирования газовых течений с учетом тепловых процессов, внешних гравитационных сил и самогравитации;
• адаптация численных алгоритмов для вычислений на компьютерах с массивно-парал-лельной архитектурой;
• исследование развития гофрировочной неустойчивости ударных волн в моделях галактического газового диска с учетом спирального потенциала звезд;
• моделирование кольцевых структур в дисковых галактиках;
• изучение возможности формирования облачной структуры в модели газовой компоненты Млечного Пути;
• исследование эффективности параллельного алгоритма на суперкомпьютере СКИФ МГУ «ЧЕБЫШЕВ»;
• сравнение результатов моделирования с наблюдениями.
Характерной мелкомасштабной особенностью большинства галактик с глобальным правильным спиральным узором являются шпуры (spurs), отходящие от спирального рукава почти перпендикулярно или пересекающие его [6]. Типичная длина этих образований лежит в пределах 100-1000 пк (рис. 1). Такие структуры принято называть feathers в смысле «оперение», «выступ», «гребень» или spurs («шпора», «отросток»), Обычно под шпурами (spurs) понимают более мощные и развитые структуры, непосредственно примыкающие к основной части спирального рукава (положению ударной волны), под оперением (feathers) - более слабые и протяженные структуры, отчасти являющиеся продолжением шпуров. Наиболее важным наблюдаемым проявлением шпуров является повышенная плотность газа и интенсивные процессы звездообразования [6].
Большинство ближайших галактик с правильной спиральной структурой демонстрирует наличие развитой системы шпуров, например, NGC 628, NGC 1232, NGC 3031, NGC 3184, NGC 4321, NGC 5194, NGC 5236, NGC 5457, NGC 4725, NGC 7424, IC0342. Изучение шпуров существенно осложняется сильной неоднородностью пыли в области спиральных рукавов,
которая одновременно выступает в качестве одного из индикаторов самих шпур. Но в целом нарушение гладкости спирального узора следует считать характерной особенностью спиральных галактик, хотя количественные различия существенны, и трудно ожидать, что все они могут быть объяснены одним универсальным физическим механизмом.
МОДЕЛЬ И ЧИСЛЕННЫЙ МЕТОД
Динамика газового галактического диска в общем случае описывается трехмерной системой уравнений гидродинамики [1]
^ + V(pM) = 0, (1)
ot
+ У(ри -и) = -Vp - pVY, (2)
ot
~\ 2 2
пи и
f (— + е) + V(p(— + £)) = е+ + Q-, (3)
at 1 I
где s = р / р(у - 1) — удельная энергия газа, р -объемная плотность газа, р - давление, и - вектор скорости, Q+, О — источники, описывающие некоторые процессы нагрева и охлаждения в газе соответственно.
Для самогравитирующего вещества с плотностью р гравитационное поле определяется уравнением Пуассона
АЧ> = 4jiGp. (4)
Газовый диск находится во внешнем гравитационном потенциале звездного диска, масса которого примерно в 10 раз больше массы газа. Учитывалась также сфероидальное темное гало с массой в 1,5-3 раза превышающее массу звездного компонента [8].
Для решения системы трехмерных уравнений газодинамики с учетом самогравитации, нагрева и охлаждения была реализована явная трехмерная TVD (Total Variation Diminishing) схема типа MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) [9]. Этот подход удовлетворяет условию невозрастания полной вариации, сохраняет монотонность численного решения для рассматриваемой задачи, характеризуется отсутствием нефизич-ных осцилляций на разрывных решениях [10], отсутствует необходимость вводить искусственную вязкость. Численный метод является эйлеровым и реализован нами для декартовой и цилиндрической систем координат. Использовался метод расщепления по пространственным переменным [11].
б
1351 ‘
Є
Рис. 1. а - оптическое изображение галактики М51; б - изображение М51 с выделением шпуров; в - отдельные участки М51 со шпурами и оперением [5]; г - пример результата компьютерного моделирования шпуров в галактических дисках
Для пространственной дискретизации был выбран не конечно-разностный подход, а более эффективный интегро-интерполяционный (метод конечных объемов), в этом случае законы сохранения выполняются на сеточном уровне. При этом значения сеточных функций заменяются средними значениями по объему ячеек, а производные вычисляются по функциям на границах ячеек. Конечно-объемная аппроксимация физических величин также позволяет избегать особенностей в центре цилиндрической системы координат при г = 0. В работе используется третий порядок аппроксимации величин по пространственным координатам. Интегрирование по времени целесообразно вести с помощью быстрого явного метода, который не понижал бы точность численной схемы. Был выбран метод Рунге-Кутта второго порядка точности.
Для вычисления потоков физических величин через границы ячеек будем использовать решение задачи Римана, на которой основаны так называемые методы Годунова [10]. Использование точного решения задачи Римана является крайне ресурсоемким, поскольку требует в каждой ячейке расчетной сетки решения нелинейного алгебраического уравнения. Воспользуемся приближенным подходом, основанным на методе Хартена - Лакса - ван Лиира (ИЬЬС), который позволяет одновременное учитывать наличие ударных волны, контактных и тангенциальных разрывов [12]. Данный метод был модифицирован для сквозного расчета границы «газ-вакуум» [13], что является необходимым для моделирования астрофизических систем, в которых могут возникать области крайне низкой концентрации вещества.
I
Рис. 2. Демонстрация чувствительности
результатов расчетов к пространственному разрешению: сравнение результатов
расчетов на декартовой координатной сетке с различным пространственным
разрешением (256^256 и 2048x2048). Аналогичный эффект имеем для цилиндрической геометрии
Получение реалистичных результатов численных расчетов требует высокого пространственного разрешения (рис. 2).
Необходимость увеличения числа ячеек для газодинамической модели приводит с одной стороны к уменьшению шага интегрирования по времени, с другой - к росту количества операций для расчета большего числа узлов сетки. Только использование параллельных вычислений на суперкомпьютерах дает возможность моделировать с высоким пространственным и временным разрешением [14].
РАСПАРАЛЛЕЛИВАНИЕ ЧИСЛЕННОГО АЛГОРИТМА
Распараллеливание численного кода выполнялось с помощью совместного использования стандартов МР1 и ОрепМР [15]. Заданная расчетная область распределялась между процессорами [16]. Декомпозиция на подобласти проводилось вдоль всех трех координатных осей: вдоль радиуса, вдоль угла и вдоль вертикально-
а
г
го направления в случае цилиндрической системы координат (см. рис. 3).
Для получения решения во всей расчетной области необходимо обеспечивать непрерывную связь соседних процессов за счет обмена данными в фиктивных граничных ячейках. Поскольку разработанная численная ТУБ схема имеет третий порядок аппроксимации по пространству, то соседним процессам необходимо обмениваться всеми физическими величинами, находящимися в двух приграничных слоях ячеек. В силу того, что данные для граничных ячеек расположены в памяти непоследовательно, то вдоль радиуса можно пересылать границу целиком (плоскость), вдоль угла последовательность одномерных массивов и вдоль вертикальной координаты - поэлементно. Таким образом, время обмена границами на каждом гидродинамическом шаге вдоль координатных осей относится как 1:2:10.
Для того, чтобы компенсировать потерю эффективности при пересылке граничных условий вдоль вертикальной координаты, по этому направлению декомпозиция должна быть наименее подробной.
Рис. 2. Декомпозиция расчетной области вдоль трех направлений для цилиндрической системы координат
Дальнейшей модификацией обмена границами стало последовательное объединение данных, необходимых для отправки, и их дальнейшая пересылка. Для вычислительных кластеров, состоящих из многоядерных процессоров, эффективным оказалось использование стандарта ОрепМР внутри каждого процесса. В этом случае происходит ускорение из-за отсутствия необходимости производить обмен фиктивными ячейками.
Структура численного газодинамического кода позволяет распараллелить вычисления с помощью ОрепМР в течение всего времени расчета на каждом процессе. В численном гид-
родинамическом алгоритме также проведена оптимизация распараллеливания за счет использования неблокирующей пересылки данных стандарта МР1 [15].
Алгоритм расчета самогравитации является более ресурсоемким по сравнению с газодинамикой, поскольку на каждом шаге интегрирования требуется адаптировать решение уравнения Пуассона под параллельные вычисления. Кратко опишем алгоритм нахождения гравитационного взаимодействия в газе TreeCode [17]. Совокупность узлов расчетной сетки может быть организована в иерархическую систему групп в форме древовидной структуры. Для сортировки частиц используется иерархия кубов, каждому узлу соответствует 8 потомков. В самой верхней части дерева находятся все частицы. Деление кубов продолжается до тех пор, пока в кубе не останется один узел, после чего для каждой ячейки или узла в дереве необходимо найти общую массу и центр масс. Для более точного вычисления гравитации была реализована возможность учитывать квадрупольные слагаемые разложения гравитационного поля.
Эффективность распараллеливания численного алгоритма была оценена в серии экспериментов (с числом ядер п = 8, 10, 16, 32, 64, 80, 100, 160, 200, 400, 800) для двух моделей: с учетом термодинамических процессов и без них (рис. 5, 6). Из графика следует, что наибольшая эффективность (около 70 %) достигается при одновременном использовании в расчетах не более 200 ядер. Графики ускорения построены для тех же наборов экспериментов (рис. 5, 6). Ускорение значительно выше в моделях, учитывающих тепловые процессы, чем для адиабатических моделей. Это обусловлено существенно большим объемом дополнительных вычислений, не влияющих на количество пересылаемых данных.
Вычисления проводились на суперкомпьютере СКИФ МГУ «ЧЕБЫШЕВ» (производитель Т-Платформы). Краткая характеристика кластера: пиковая производительность 60 ТИор^, 1250 процессоров, 5000 ядер в системе, объём оперативной памяти 5,5 Тбайт, модель процессора Щ;е1 Хеоп Е5472 3,0 ГГц.
Далее рассмотрим основные результаты численного моделирования, полученные с использованием созданного программного продукта.
Рис. 3. Эволюция поверхностной плотности газового диска с учетом неосесимметичного гравитационного потенциала звездного компонента. На начальных этапах происходит формирование системы ударных волн. В дальнейшем под воздействием гофрировочной неустойчивости на нелинейной стадии развивается сложная система протяженных шпуров
Рис. 4. Графики ускорения параллельных расчетов для моделей с учетом и без учета тепловых процессов
Рис. 5. Графики эффективности параллелизма для реализованного численного алгоритма с учетом и без учета тепловых процессов в газе
РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ
Гофрировочная неустойчивость в газовом диске
Рассмотрим возможность генерации шпуров за счет чисто гидродинамических эффектов, связанных с существенно нелинейной стадией развития гофрировочной (wiggle) неустойчивости фронта глобальной ударной волны в газовом диске при наличии спиральной волны плотности звездного компонента [18].
В модели имеется 8 свободных физических параметров, варьирование которых в широких пределах позволило рассмотреть формирование различных структур, связанных с развитием гофрировочной неустойчивости глобальной ударной волны. Значительное количество свободных параметров задачи приводит к необходимости проведения большого числа компьютерных экспериментов для исследования одной модели галактического газового диска. Для одной модели проводилось до 10-15 расчетов
с различными наборами параметров, что занимало в среднем 10000 процессор/часов.
Возможность возникновения и скорость нарастания шпуров существенно зависят от их положения на спиральном рукаве. Как правило, очаг первоначального формирования нелинейной стадии гофрировки локализован на небольшом участке ударной волны (рис. 4). По мере увеличения длины отростка газовой спирали он может достигать соседнего участка витка спирали, индуцируя возмущения на нем. Такой механизм приводит к формированию наиболее интенсивных и протяженных шпуров с мощным оперением.
Моделирование облачной структуры Млечного Пути
Наблюдательные данные о газопылевых комплексах (ГПК) нашей Галактики показывают, что эти объекты концентрируются к спиральным рукавам Галактики, при этом свойства комплексов зависят от положения в Галактике [19].
Численные модели, описывающие динамику газа в Галактике, должны адекватно описывать мелкомасштабную структуру и внутреннюю кинематику молекулярных облаков в галактическом диске. Это предъявляет особые требования к пространственному и временному разрешению, а также адекватный учет наиболее важных физических факторов, связанных с нагревом и охлаждением.
Основываясь на результатах численных экспериментов, проведенных на суперкомпьютере СКИФ МГУ «ЧЕБЫШЕВ», удалось продемонстрировать сценарии образования облачной структуры газа в окрестностях спиральных рукавов Галактики. Достигнутое пространственное разрешение позволило изучать отдельные мелкомасштабные газовые комплексы.
В моделях, учитывающих тепловые процессы (нагрев - охлаждение), хорошо выделяются газопылевые комплексы. В их плотных частях температура составляет 50-200 К. При этом дисперсии турбулентных скоростей для таких образований порядка 3-5 км/с. Параметры этих структур, полученные в численных экспериментах, соответствуют данным наблюдений гигантских молекулярных облаках и их окружении.
В рамках численных гидродинамических расчетов продемонстрирована возможность формирования ГПК, содержащих структуры типа гигантских молекулярных облаков. Структура и кинематика газа в построенных моделях
качественно отражает наблюдаемые особенности распределения газопылевых комплексов.
Образование кольцеобразных галактических структур
Известно значительное число галактик, у которых на некотором расстоянии от центра наблюдается кольцеобразная спиралевидная структура.
Рис. 6. Слева: кольцевая структура в численном эксперименте. Справа: пример галактики N00 1512 с кольцом
Можно формально выделить три типа таких колец, не затрагивая различные механизмы их образования: ядерные кольца (галактики N00 1512, N00 4314, N00 6782), кольца на концах бара (N00 1512, N00 3081, N00 4725, N00 5905) и периферийные кольца в дисковой компоненте (М31, N00 7742, АМ 0644-741).
В данной работе рассмотрен гидродинамический механизм формирования кольцеобразных структур в галактических дисках.
На рис. 7 показан пример расчета динамики газового диска с образованием кольцеобразных структур во внешней области системы, причем во внутренней части диска имеется правильный двухрукавный узор со слабыми проявлениями нерегулярности. Напротив, кольцо обнаруживает все признаки наличия шпуров и оперения различной степени интенсивности.
ЗАКЛЮЧЕНИЕ
В рамках данной работы было проведено около 200 численных экспериментов на вычислительном кластере СКИФ МГУ «ЧЕБЫШЕВ». На эти расчеты было затрачено порядка 200000 процессор/часов. Благодаря использованию ресурсоемких вычислений на суперкомпьютере были промоделированы тонкие структуры в газовых дисках галактик. Результаты теоретических расчетов были сопоставлены с реально наблюдаемыми структурами в спиральных галак-
тиках. Сформулируем наиболее важные результаты:
• Создан универсальный программный продукт для моделирования газодинамических процессов различной геометрии с учетом само-гравитации, тепловых процессов и внешних гравитационных полей. Данный пакет программ адаптирован для использования на ЭВМ с массивно-параллельной архитектурой.
• Численные эксперименты позволили объяснить ряд наблюдаемых структур (шпуры, оперение) которые являются результатом развития гофрировочной неустойчивости ударных волн в газовом диске.
• Предложен и численно-экспериментально изучен гидродинамический механизм образования галактических колец в газовой компоненте галактик.
• Продемонстрирован процесс формирования облаков в нашей Галактике в результате развития газодинамической и тепловой неустойчивостей.
СПИСОК ЛИТЕРАТУРЫ
1. Фридман А. М., Хоперсков А. В. Физика галактических дисков. М.: Физматлит, 2011. 640 с.
2. Бочкарев Н. Г. Основы физики межзвездной среды. М.: Изд-во МГУ, 1991. 352 с.
3. Засов А. В., Сильченко О. К. Диски галактик и их эволюция // УФН, 2010. Т. 180. С. 434-439.
4. High resolution simulations of unstable modes in a collisionless disc / A. V. Khoperskov [et al.] // Astronomy & Astrophysics. 2007. V. 473. P. 31-40.
5. Detection of Dense Molecular Gas in Interarm Spurs in M51 / S. Corder [et al.] // Astrophysical Journal. 2008. V. 689. P. 148-152.
6. ASTE CO (3-2) Mapping Toward the Whole Optical Disk of M 83: Properties of Inter-arm Giant Mo-lecular-Cloud Associations / K. Muraoka [et al.] // Astrophysical Journal. 2009. V. 706. P. 1213-1225.
7. Roberts W., Hausman M. Spiral structure and star formation. I - Formation mechanisms and mean free paths // Astrophysical Journal. 1984. V. 277. P. 744-767.
8. Numerical modeling of the vertical structure and dark halo parameters in disc galaxies / A. V. Khoperskov [et al.] // Astronomische Nachrichten. 2010. V. 331. P. 731-745.
9. van Leer B. A second-order segcul to Godunov's method // Journal of Computational Physics. 1979. V. 32. P. 101-136.
10. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физмалит, 2001. 608 c.
11. Воеводин В. В., Воеводин Вл. В. Параллельные вычисления. СПб.: БХВ-Петербург, 2002. 608 с.
12. Shengtai Li. An HLLC Riemann solver for magneto-hydrodynamics // Journal of Computational Physics. 2004. V. 203. P. 344-352.
13. Еремин М. А., Хоперсков А. В., Хоперсков С. А. Конечно-объемная схема интегрирования уравнений гидродинамики // Известия ВолГТУ. Актуальные проблемы управления, вычислительной техники и информатики. 2010. Т. 13. C. 24-27.
14. Васильев В. А., Ницкий А. Ю. Исследование масштабируемости задач вычислительной гидроаэродинамики на различных многоядерных и многопроцессорных архитектурах // Вестник УГАТУ. 2010. Т. 14. C. 126-132.
15. Антонов А. С. Введение в параллельные вычисления, М.: Изд-во МГУ, 2002. 69 c.
16. Кузнецов П. В., Кайгородов О. А. Адаптация схемы Ошера для параллельных вычислений. М.: ИПМ им. М. В. Келдыша. 2002. Т. 59. 10 с.
17. Barnes J., Hut P. A hierarchical O (N log N) force-calculation algorithm // Nature. 1986. V. 324. P. 446-449.
18. Wada K. Instabilities of spiral shocks - II. A quasi-steady state in the multi-phase inhomogeneous ISM // Astrophysical Journal. 2008. V. 675. P. 188-193.
19. Kalberla P. M. W., Kerp J. The Hi Distribution of the Milky Way // Annual Review of Astronomy & Astrophysics. 2009. V. 47. P. 27-61.
ОБ АВТОРАХ
Хоперсков Сергей Александрович, магистрант каф. теоретическ. физики и волновых процессов Волгоградск. гос. ун-та.
Хоперсков Александр Валентинович, проф., зав. каф. инф. систем и компьютерн. моделирования, Волгоградск. гос. ун-та. Д-р физ.-мат. наук.
Засов Анатолий Владимирович, проф., руководитель отдела Внегалактической астрономии гос. ас-трономическ. ин-та им. П. К. Штернберга МГУ. Д-р физ.-мат. наук.
Бутенко Мария Анатольевна, асп. каф. инф. систем и компьютерн. моделирования, Волгоградск. гос. ун-та.
