CC BY
24
ISSN 1992-6502 ( P ri nt)_
2014. Т. 18, № 1 (62). С. 230-240
Ъыьмт QjrAQnQj
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 519.6
Параллельный программный комплекс
для моделирования астрофизических струйных выбросов
на гибридных вычислительных системах с общей памятью
12 3
м. п. Галанин , в. в. Лукин , в. м. Чечеткин
1 galan@keldysh.ru, 2 vvlukin@gmail.com
1-3 ФГБУН «Институт прикладной математики им. М. В. Келдыша РАН» 1,2 ФГБОУ ВПО «Московский государственный технический университет им. Н. Э. Баумана»
Поступила в редакцию 22.12.2013
Аннотация. Рассмотрена модель образования, коллимации и радиационного ускорения плазменного струйного выброса в окрестности компактного объекта. Модель включает в себя систему уравнений радиационной магнитной гидродинамики в двумерном цилиндрически симметричном приближении. Разработан параллельный программный комплекс для систем с общей памятью, использующий технологии программирования OpenMP и nVidia CUDA, реализующий разностную схему на треугольных неструктурированных сетках для решения уравнений МГД и метод дискретных направлений решения уравнения переноса излучения. Показано, что ускоряемое радиационным давлением вещество достигает высоких скоростей, соответствующих наблюдаемым в астрофизических системах.
Ключевые слова: радиационная магнитная гидродинамика; вычисления на графических процессорах; астрофизические джеты; системы с общей памятью.
Одним из наиболее интересных классов астрофизических процессов является образование струйных выбросов, джетов (от англ. jet -струя), формирующихся в ядрах активных галактик (например, эллиптической галактике M87), микроквазарах (например, двойной звездной системе SS433 [1, 2]), протозвездных объектах и ряде других систем. Тонкие (поперечный размер обычно не превышает нескольких парсек) биполярные струи вещества около таких объектов тянутся на сотни и тысячи световых лет, заканчиваясь гигантскими облаками газа.
На настоящий момент известно достаточно большое количество джетов, большинство из них наблюдается только в радиодиапазоне. Размер джета, истекающего из ядра галактики M87, достигает 5000 световых лет. Струя
Работа выполнена при финансовой поддержке РФФИ (проекты 12-01-00109, 14-01-31496, 12-0212096, 12-02-00687), гранта Президента Российской Федерации для государственной поддержки ведущих научных школ Российской Федерации (проект НШ-6061.2014.2). Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2012».
состоит из быстро движущихся заряженных частиц, сконцентрированных в узлы размером до 10 световых лет (см. рис. 1), и имеет вид конуса с углом раствора около 6°. Скорость течения вещества в джете галактики M87 достигает 0.8с, где с - скорость света, скорость вещества в джете SS433 - 0.26с .
4 July 2003
Рис. 1. Результаты наблюдений объекта SS433 на телескопе УЬВА [2]: смена максимумов и минимумов яркости в указанных точках
____г
Рис. 2. Схема модели системы, порождающей джет, и расчетная область
В данной работе построена и исследована математическая модель ускорения вещества джета в канале над горячим гравитирующим объектом под действием давления излучения тонкого диска, окружающего компактный объект. Модель включает в себя систему уравнений идеальной радиационной магнитной гидродинамики (МГД) в двумерном цилиндрически симметричном приближении и основана на МГД-модели образования ускоряюшего канала [3].
Моделирование излучения с помощью уравнения переноса с учетом интеграла рассеяния представляет вычислительно тяжелую задачу. В двумерном случае для адекватного расчета поля излучения применены параллельные алгоритмы для высокопроизводительных вычислительных систем с общей памятью, включающих графические ускорители.
ПОСТАНОВКА ЗАДАЧИ
Предлагается использовать различные эффекты для объяснения различных наблюдаемых свойств джетов. Так, для объяснения высокой степени коллимации джета принято использовать модели магнитной гидродинамики (МГД), в которых коллимация потока
2.5 2
1.5
г
1
0.5 О
происходит под действием осевого и частично азимутального компонентов магнитного поля
[3, 4].
В работе [3] показано, что над гравити-рующим объектом, окруженным тонким вращающимся диском и погруженным в облако аккрецирующей плазмы (см. рис. 2), образуется замагниченный канал, внутри которого развивается струйный выброс плазмы, источником которой является тонкий диск (см. рис. 3).
Полученный канализированный выброс оказывается устойчивым во времени, канал имеет плотные, оптически толстые стенки, в то время как вещество внутри джета сильно разрежено, его оптическая толща (для случая томпсоновского рассеяния) составляет
т = атпЦ «6.7-10Л (1)
где Х0 - пространственный масштаб задачи, п -концентрация вещества в канале, ат -томпсоновское сечение рассеяния.
В работе [5] рассмотрена нульмерная динамическая модель ускорения твердотельного сгустка в пустотелом канале давлением излучения горячего «дна» канала, моделирующего окрестности компактного объекта. Показано, что в подобной модели сгусток достигает субсветовых скоростей, вплоть до 0.9с .
Рис. 3. Ускоряющий канал в модели [5]: распределения плотности и скорости в момент времени t = 15
(установившийся режим выброса)
Мы будем использовать схему модели, изображенную на рис. 2, и результаты работы [3] в качестве начальных условий. Рассмотрим квазистационарную моноэнергетическую модель распространения излучения и модель магнитной гидродинамики с идеальной проводимостью плазмы в двумерном цилиндрически симметричном приближении. Пусть вещество системы представляет собой невязкий совершенный идеально ионизованный газ, электрическое сопротивление среды отсутствует. Магнитное поле вморожено в тонкий идеально проводящий диск и вращается вместе с ним, приобретая коллимирующий осевой и ускоряющий азимутальный компоненты. Гравитационное поле определяется гравитацией центрального тела (звезды), самогравитация газа не учитывается. Фотоны испытывают однократное рассеяние на электронах (томп-соновское рассеяние), источником излучения в модели служит только тонкий диск, излучение сфокусировано в ускоряющий канал.
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Система уравнений. Полная система уравнений радиационной магнитной гидродинамки в рассматриваемом приближении имеет вид [6]:
Яр
Ы
+ Уру = 0,
|(ру + С ) + У-(п + Т)=
= — (Ухб)хВ + Е , 4лV ' '
я
—(в + и)+У-(у(в + р) + W) =
Я
= — ((УхВ)хВ)-у + е • у,
4л
(2)
(3)
(4)
ш • У1 х, га) + к (¿, х)1 (¿, х, га) =
= Р(^, х) |Г(^ х, ш, ш')1 х, ш')^ш', (5)
— = У х (у х В),
(6)
где р - плотность плазмы, у - вектор скорости
вещества, II - тензор плотности потока импульса вещества, П = рд + ру^, 5,у -
символ Кронекера, р - газовое давление, в -плотность энергии вещества, В - вектор магнитной индукции, Е - объемная плотность
гравитационной силы, I - интенсивность излучения, Г(^ х, ю • ю') - индикатриса рассеяния, к(£, х) - коэффициент ослабления,
к = а + р, а(£, х) - коэффициент поглощения, р(£, х) - коэффициент рассеяния излучения в
веществе, С = W/c2 - плотность имульса
излучения, W = \шЫш - поток энергии
и
излучения, и = — \Ыш - плотность энергии
С и
излучения, = — ^¡щЫш - тензор плотности
С и
потока импульса излучения. Газ совершенный и имеет уравнение состояния р = рв(у —1), где в - удельная внутренняя энергия вещества, откуда
в =
_р| у
- + -
2 у — 1
Основную трудность при исследовании этой модели представляет уравнение (5) - уравнение переноса излучения (УПИ). Это связано с принципиальной многомерностью уравнения: вообще говоря, интенсивность I зависит от трех пространственных координат (радиус-вектор х), двух угловых координат (единичный вектор ю можно задать с помощью двух координат ^ = ео80 и ф, где 0 и ф -соответственно полярный и азимутальный углы) и времени t. Последняя зависимость является в квазистационарной модели переноса излучения неявной и выражается только в зависимости радиационных коэффициентов вещества от термодинамических параметров плазмы.
Отметим, что влияние поля излучения на движение вещества описывается в уравнениях (3) и (4) с помощью градиентов различных моментов интенсивности излучения. Данное обстоятельство при численном решении системы (2)-(6) приводит к необходимости получать насколько возможно гладкие распределения интенсивности, минимизируя эффекты численной немонотонности, либо получать сразу указанные градиенты моментов. Подобное обстоятельство становится критическим, например, для методов класса Монте-Карло, поскольку эти методы обладают высокой статистической дисперсией.
В качестве индикатрисы рассеяния будем использовать рэлеевскую индикатрису
4
а сечение рассеяния ит = 6.652х 10—29 см2 (коэффициент рассеяния, соответственно, равен р = пат, где п - концентрация вещества).
Мы предполагаем, что в начале координат находится тело массы М, являющееся источником гравитационного поля. Чтобы
Г^, х, ш, ш') =—(1 + (ш • ш')2),
2
и
избежать сингулярности в начале координат, будем задавать гравитационную силу
следующим образом (Я = л/г2 + г2 ):
р = -о Мр г,
г Я2 Я
Т7 ПМР
= -о—г
г
р = -оМр г
Я1 Я
Я > г,
р = -О МР г, Я < г*, г
где О - гравитационная постоянная.
Система уравнений стандартным образом приводится к безразмерному виду, в качестве основных масштабов выбираются масштабы расстояния Х0, плотности р и времени ^.
Безразмерный параметр g = GMt2/L>0 =0.5.
Граничные условия. Итак, предполагаем, что в пространстве, заполненном идеально проводящей плазмой, имеются диск и гравитирующее тело в центре. Расчетная область [0,гм ] х [0, гм ] представлена на рис. 2.
На внешней цилиндрической границе области будем задавать условие сверхзвукового сферического втекания незамагниченной незакрученной межзвездной плазмы.
Верхняя граница расчетной области моделирует переход потока к режиму течения на бесконечности [7].
На оси вращения системы г = 0 ставится условие ограниченности решения.
Нижняя граница разбивается на две части. На части границы г = 0, г < г < ^, где ^ -радиус тонкого диска, ставятся условия, соответствующие экваториальной симметрии, а на части г = 0 , 0< г < г моделируется тонкий диск. Диск вращается со скоростью ш(г) = ш- (1- (г/г )2) . В вещество диска вморожено магнитное поле, имеющее только осевой компонент Вг0 (г). Диск считается идеально проводящим. Он является источником вещества джета - вещество поступает в расчетную область со скоростями, определяемыми течением над диском.
Кроме того, необходимо поставить граничные условия для поля излучения, смоделировав источник в окрестности центрального грави-тирующего объекта. Будем предполагать, что эти окрестности излучают с интенсивностью абсолютно черного тела, имеющего характерную температуру 7 -104 К, что отвечает наблюдательным данным об объекте SS433 [1]. При этом излучение сфокусировано внутрь ускоряющего канала: излучающей является граница г = 0,0ПгЛ 0.2, излучение распро-
страняется вдоль направлений, для которых cos9 > 0.9, где 0 - полярный угол луча.
ЧИСЛЕННЫЙ МЕТОД И ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
Конечно-разностный алгоритм. Для
решения системы уравнений идеальной радиационной магнитной гидродинамики в двухмерной цилиндрически симметричной постановке использованы треугольные неструктурированные сетки и разностный метод, основанный на расщеплении по физическим процессам. Пересчет неизвестных величин в разностных ячейках на каждом временном шаге состоит из четырех этапов.
1. Решение газодинамической системы уравнений методом типа Годунова (HLLC), а также учет геометрии задачи и действия гравитационного поля в виде слагаемых в правой части [8].
2. Аппроксимация уравнения Фарадея на разностной ячейке в цилиндрической системе координат, согласованная с геометрией задачи.
3. Интегрирование УПИ методом дискретных направлений (МДН) с дискретизацией сферы направлений «от границы» [9].
4. Восполнение газовых переменных за счет действия магнитных и радиационных сил.
Реализация данного алгоритма на треугольных неструктурированных сетках в целом изложена в [9].
Особенности реализации МДН в случае цилиндрической симметрии. Как уже отмечалось, с вычислительной точки зрения наиболее требовательной к ресурсам является процедура решения уравнения переноса излучения (5). Для решения этого уравнения выбран метод дискретных направлений. В соответствии с данным методом необходимо провести интегирование УПИ вдоль каждого из выбранных дискретных направлений, сведя тем самым многомерное интегродифференциальное уравнение (5) к множеству одномерных обыкновенных дифференциальных уравнений (в случае если интеграл рассеяния учитывается итерационно):
<1 (ц, га)
<ц
■ + к (ц)1 (Ц, га) = £ (ц, га),
где ц - параметр, изменяющийся вдоль луча,
проходящего через данный узел пространственной разностной сетки, и имеющий смысл расстояния, £(ц, га) - источник рассеянного излучения.
Таким образом, для каждой точки пространственной сетки х; необходимо про-
вести трассировку лучей, приходящих в нее вдоль векторов га. выбранной с учетом удаленности точки х; от границы сферы направлений.
Отметим, что аккуратный учет граничных условий для поля излучения приводит к выбору своего набора дискретных направлений для каждой пространственной точки, что серьезно увеличивает требования к вычислительным ресурсам, прежде всего к памяти.
Отметим также, что в отличие от чисто двумерной задачи, наличие цилиндрической симметрии приводит к необходимости трассировать направления распространения излучения в трехмерной области (рис. 4). При этом в силу симметрии достаточно вычислить интенсивность излучения I(х, га) = (г, г, га) в полуплоскости гОг.
Для вычисления интенсивности /г(> (г, г, га) в соответствии с методом дискретных направлений необходимо протрассировать луч, проходящий
через точку [г,0,г]т в направлении га, по полуцилиндрической трехмерной области А , причем треугольные ячейки плоской разностной сетки в области О (залита серым цветом на рис. 3 слева) соответствуют пространственным кольцам треугольного сечения, которые заполняют А . Таким образом, получаем задачу трассировки прямолинейного луча в области А , состоящей из кольцевых подобластей с постоянными оптическими свойствами вещества (пример кольцевой подобласти изображен на рис. 3, справа). Эту задачу легко свести к трассировке кривой второго порядка (гиперболы) в плоской области О.
Отметим, что трассировка гиперболы на треугольной сетке является существенно более затратной процедурой (затраты машинного
времени возрастают на порядок), чем трассировка прямой, поэтому предпочтительным является хранение всей информации о трассировке лучей в оперативной памяти в процессе динамического расчета с интегрированием квазистационарного УПИ на каждом временном слое. Подобный объем памяти на данный момент не доступен персональным компьютерам, но вполне достижим на суперкомпьютерной вычислительной технике.
Учет рассеяного излучения. Учет интеграла рассеяния
£ (г, х, га) = |Г(г, х, га, га')/(г, х, га' )<га' (7)
и
в примененной численной схеме производится итерационно в соответствии с известными схемами для моделирования однократного рассеяния. На каждом шаге по времени для расчета интенсивности излучения производится две итерации. На первой из них решается уравнение переноса излучения в виде
га - VI(1) (г, х, га) + к (г, х)/(1) (г, х, га) = 0,
т. е. в предположении нулевого вклада рассеянных фотонов в общую интенсивность.
Далее вычисляется функция источника рас-сеяного излучения
£(1) (г, х, га) = /Г(г, х, га, га') I(1) (г, х, га')<га'
и
и производится пересчет поля излучения в расчетной области путем интегрирования уравнения
га - VI(2) (г, х, га) + к (г, х)/(2) (г, х, га) =
= р(г, х)£ (1)(г, х, га)
с известной правой частью.
Подобный подход требует на каждом шаге
Рис. 4. Цилиндрическая симметрия и трассировка луча в задаче о переносе излучения
по времени в каждом узле пространственной сетки производить интегрирование функции интенсивности I(1)(t, x, w) по всей сфере направлений. Взятие таких интегралов является наиболее затратной в вычислительном плане процедурой, причем тем более тяжелой, чем более подробно разбиваются сферы направлений в узлах сетки. В то же время такие интегралы берутся независимо в разных пространственных точках, алгоритм их взятия не содержит управляющих конструкций ветвления, что позволяет широко использовать параллельные вычислительные технологии, в частности технологию CUDA для вычислений на графических ускорителях.
Параллельные вычислительные технологии. Для решения ресурсоемких задач, подобных рассматриваемым в данной работе моделям радиационной МГД, в настоящее время широко применяется суперкомпьютерная техника, реализующая параллельные вычислительные технологии для систем с общей и распределенной памятью [10]. Мы уже указали на высокие требования метода дискретных направлений в целом и процедуры трассировки лучей и интегрирования вдоль них УПИ в частности к ресурсам памяти вычислительного модуля. В данном случае наиболее целесообразным оказывается именно использование SMP-моду-лей с максимальным доступным на данном узле кластера числом процессорных ядер и максимальной доступной памятью. В случае кластера К-100 ИПМ им. М. В. Келдыша РАН имеется возможность использования узла кластера целиком в качестве одного вычислительного модуля, содержащего 12 процессорных ядер, 96 Гб оперативной памяти и 3 графических ускорителя nVidia Tesla.
Численный метод реализован в виде программного комплекса на языках программирования Фортран-90 и C++. Программный комплекс разработан для применения на высокопроизводительных кластерах с графическими ускорителями и реализует следующий подход к распараллеливанию на основе технологий для систем с общей памятью.
Применение технологии OpenMP. Интегрирование УПИ вдоль протрассированных и сохраненных в ОЗУ лучей, соответствующих
данному узлу пространственной сетки, осуществляется независимо для каждого узла. При этом информация о трассе луча (номер проходимой ячейки, длина отрезка луча) используется только однажды. В данном случае удобно производить параллельные вычисления над общей оперативной памятью, используя возможности многоядерных процессоров общего назначения. В этом случае узел кластера используется как SMP вычислительный модуль с 12 вычислительными ядрами (каждое из которых включает одно процессорное ядро) и 96 Гб общей для всех ядер памяти.
Применение технологии nVidia CUDA. Будучи более тяжелой операцией, вычисление интеграла рассеяния одновременно является более удобной для распараллеливания. В рамках данной операции для каждого узла пространственной сетки вычисляется Nxat интегралов вида (7), где N^at - количество дискретных направлений, в которые производится рассеяние излучения центральной аккреционной системы модели. При этом для вычисления каждого из этих Nxat интегралов используются одни и те же данные об интенсивности приходящего излучения и телесном угле, откуда это излучение приходит в рассматриваемый пространственный узел. Именно в таких задачах наиболее эффективными оказываются вычисления на графических ускорителях благодаря относительно низкой доле времени обменов между памятью ускорителя и основной памятью узла в общем времени вычислений. В этом случае узел кластера используется как модуль с одним вычислительным ядром, включающим одно из процессорных ядер, один графический ускоритель с интегрированной памятью и 96 Гб памяти общего назначения. Вычисление всей совокупности интегралов вида (7) в данном пространственном узле производится в рамках блока тредов одного графического мультипроцессора, а каждый отдельный интеграл вычисляется в рамках одного треда в указанном блоке. При этом загрузка распределения интенсивности и телесных углов для данного пространственного узла и данной сферы направлений производится однократно, обеспечивая данные для всего блока тредов.
Рис. 5. Начальное распределение интенсивности I и модуля силы давления излучения
Использование указанных подходов позволило существенно ускорить работу программного комплекса. В частности, процедура трассировки и интегрирования УПИ вдоль луча с использованием технологии OpenMP на 12-ядерном узле кластера К-100 выполняется в 10.8 раз быстрее, чем на одном процессорном ядре, а процедура вычисления источников рассеяного излучения во всех пространственных узлах выполняется на графическом процессоре nVidia Tesla в 82.3 быстрее, чем на одном процессорном ядре (все измерения производились на одном узле кластера К-100, процессор Intel Xeon X5670 2.93 ГГц).
РЕЗУЛЬТАТЫ РАСЧЕТОВ
Вычисления производились с использованием кластера К-100 ИПМ им. М. В. Келдыша РАН. Использованы треугольные сетки, состоящие из 8000 ячеек, со сгущением у диска. Проведена серия вычислительных экспериментов для разных вариантов граничных условий на тонком диске. Подробно рассмотрим случай, когда плотность ветра с тонкого диска ограничена снизу величиной р^. Шаг по времени в расчетах выбирался автоматически в зависимости от скорости течения плазмы в расчетной области. Во избежание развития неустойчивостей число Куранта принималось равным СЕЬ = 0.05.
Рис. 6. Распределения параметров течения вещества вдоль оси г = 0: плотности р , давления р, объемной силы радиационного давления, модуля скорости |у|
Ускорение выброса в канале. Начальными условиями для расчетов служат распределения величин в момент времени г = 18.05, полученные в [3]. В рассматриваемой области существует замагниченный канал, внутри которого движется разреженное вещество струйного выброса, причем плотность вещества падает с удалением от центральной области (рис. 6).
«Включение» поля излучения приводит к возникновению дополнительной ускоряющей силы, действующей на разреженное вещество, величина этой силы тем больше, чем ближе к источнику излучения. На рис. 5 приведено начальное распределение интенсивности излучения и модуля силы радиационного давления (под силой в данном случае понимается член V- Т, поскольку член дО/дг в (3) мал по модулю и на характер течения не влияет), действующей на вещество выброса. Действие этой силы в обоих расчетах приводит к тому, что скорость течения выброса быстро достигает значений порядка 103, а ускорение происходит в ближайшей окрестности излучающего центрального объекта.
Сила радиационного давления убывает с удалением от источника излучения, так что нижние (относительно оси Ог ) слои ускоряются быстрее верхних, что приводит к сжатию и повышению давления в потоке выброса. Наконец, к моменту времени г = 18.055 наибольшая скорость течения достигает максимума за все время расчета и дальше не увеличивается, начальные эффекты считаются пройденными, устанавливается стабильный режим выброса. На рис. 6 изображены распределения основных величин течения вещества выброса вдоль оси Ог в момент времени г = 18.075. На рис. 7 изображены распределения плотности, модуля
скорости а также полоидальных мгновенных траекторий и полоидальных силовых линий магнитного поля в расчетной.
Полученные скорости достигают 1/5сй, где с^ - безразмерная скорость света. При этом выброс остается хорошо коллимированным, магнитное поле в канале сохраняет свою структуру, вокруг наиболее высокоскоростной части выброса сохраняются мягкие стенки из закрученного магнитного поля. Давление и плотность газа, начиная с некоторой высоты г, падают при удалении от центрального объекта, поэтому выброс не будет замедляться вне расчетной области и, достигнув предельного значения скорости, будет распространяться в разреженную среду, сохраняя высокую энергетику. В то же время максимум давления приходится на окрестности плоскости г = 0.4, т. е. на верхнюю часть сужения канала. Подобное распределение градиента давления вдоль оси Ог может приводить к возникновению эффекта «запирания» вещества в окрестности центрального объекта в случае, если силы радиационного давления будет недостаточно, чтобы преодолеть этот своеобразный «барьер».
В расчетах устанавливается периодический режим выброса, при котором сохраняется высокий темп течения вещества, плазма постоянно поступает в канал с тонкого диска в экваториальной плоскости, не давая джету стать более разреженным. Отметим, что для устойчивого функционирования механизма формирования и ускорения струйного выброса необходимо постоянное поступление вещества внутрь ускоряющего канала с плотностью потока не меньше некоторого предельного значения, например р^ =0.01.
Рис. 8. Обострение волны скорости в ходе ускорения всплеска и образование фронта ударной волны. Распределение модуля скорости вдоль оси Ох в последовательные моменты времени (слева) и волна
плотности в районе фронта волны (справа)
Периодический режим выброса. Особый интерес вызывает установившийся в расчетах периодический режим выброса вещества. В подобном режиме сохраняется средняя высокая скорость истечения вещества, при этом наблюдаются периодические всплески скорости потока. Период всплесков Т =0.01 совпадает со временем, необходимым плазме, чтобы уйти от центрального объекта и преодолеть сужение ускоряющего канала, выйдя на стабильный уровень скорости.
Система переходит в колебательный режим сразу после установления скорости потока, когда первая порция плазмы, ускоренной излучением, покидает расчетную область. После ухода этой порции над тонким диском образуется область разреженности вещества. В сужении канала над излучающим центральным объектом образуется барьер из плазмы высокого давления, который частично «запирает» вещество в получившейся «каверне» . Поток массы с диска ограничен снизу путем поддержания плотности плазмы не менее значения р^. Цикл действия полученного механизма выглядит следующим образом.
1. Над тонким диском существует разреженная каверна, запертая сверху барьером давления горячего газа шлейфа предыдущего всплеска. За счет преобладания потока массы с диска над потоком массы через барьер каверна начинает заполняться плазмой, становится более непрозрачной для излучения центрального объекта.
2. Повышение плотности газа в каверне приводит к увеличению силы радиационного давления на порядок. Давление излучения тем больше, чем ближе вещество находится к
центральному объекту - вещество в каверне начинает сжиматься давлением излучения снизу, давление газа в каверне растет.
3. Давление газа в каверне превышает давление газа барьера в сужении канала, накопленное в каверне плотное вещество образует быстро движущийся вверх по оси Ох горячий пузырь. Образуется всплеск скорости, проходящий сужение канала и растущий далее за счет расширения пузыря.
4. Пузырь выходит в расширяющуюся часть канала, плотность его падает вместе с расширением канала, скорость продолжает нарастать, в распределении скорости образуется фронт всплеска, причем нижние слои газа движутся быстрее, накатывая на верхние и заостряя фронт (рис. 8). За пузырем образуется область разреженного газа низкого давления, барьер давления в сужении канала закрывается, начинается накопление новой порции выбрасываемого вещества.
5. На выходе из расчетной области фронт всплеска формирует ударную волну, нагребающую на себя газ из фонового течения в канале. Формируется сгусток вещества, движущегося со скоростями порядка 1/5С , где С - безразмерная скорость света. Сгусток продолжает расширяться и покидает расчетную область. Каверна над центральным объектом заполняется плазмой, истекающей с диска.
Обсуждение результатов. Перейдем от результатов расчетов, представленных в безразмерном виде, к оценкам значений размерных переменных задачи. Примем следующие масштабы величин: концентрация п0 = 108 см—; плотность
р0 =3.34х10~
см
(молекулярный водород); линейный размер Ь = 1015 см ; напряженность магнитного поля
B =0.06 Э;
масса
->33
центрального
тела
М = 3МС =6 х 10 г; временной масштаб задачи г0 = 1. 12х109 с; масштаб скорости
V =0.9х106 см/с.
Расчет моделирует истечение от формирующейся протозвезды с массой М = 3М , окруженной околозвездным диском радиусом г =0.6Ь0 «40а.е. Диск пронизан полои-дальным магнитным полем, напряженность которого составляет 0.06 Э. На систему аккре-цирует сверхзвуковой поток вещества с темпом аккреции 5 х 105 Мс/год. Центральные области системы излучают, интенсивность излучения соответствует излучению абсолютно черного тела температурой 7 х 104 К . Над звездой с диском сформирован замагниченный канал, содержащий разреженное вещество, источником которого является диск. Излучение центрального объекта сфокусировано внутрь канала, претерпевает лишь однократное рассеяние на электронах (томпсоновское рассеяние). Перпендикулярно экваториальной плоскости диска формируется коллимированное истечение вещества (джет), ускоряемое давлением излучения. Скорость потока в среднем составляет 2х104 км/с, угол раствора джета « 10°. Поток состоит из сгустков вещества, движущихся со скоростью большей, чем скорость фонового течения. Максимальная скорость сгустков достигает 5х104 км/с. Период выброса сгустков составляет 13 дней.
На образование сгустков принципиальное влияние оказывают три фактора: радиационное давление излучения центрального объекта, эффективно ускоряющего плотное вещество и мало действующего на разреженную плазму; образование на начальном этапе радиационного ускорения выброса в сужении канала «барьера» газового давления, приводящего к запиранию плазмы в каверне над излучающим диском; ограниченность снизу потока вещества с диска, что обеспечивает постоянное пополнение каверны над излучающим диском веществом и воспроизводство рассмотренного механизма.
Полученный струйный выброс является устойчивым при условии ограничения снизу потока вещества с диска, которое препятствует «выдуванию» плазмы из канала и прекращению действия силы радиационного давления. Кол-
лимация джета обеспечивается осевым магнитным полем во внутренней части канала и наличием у канала мягких спиральных стенок из закрученного магнитного поля.
ЗАКЛЮЧЕНИЕ
Построена и исследована модель радиационного ускорения вещества в замагниченном канале над тонким идеально проводящим диском. Математическая модель состоит из уравнения переноса излучения и системы уравнений идеальной магнитной гидродинамики с учетом гравитационной силы, создаваемой центральным объектом, и давления излучения. Для численного решения системы уравнений применена модификация алгоритмов, разработанных и рассмотренных в [3]. Рассмотренные методы (разностная схема для системы уравнений идеальной МГД и метод дискретных направлений для УПИ) адаптированы для расчетов в цилиндрической системе координат и реализованы в виде программного комплекса для высокопроизводительных SMP-машин с графическими ускорителями. Построенный программный комплекс показал высокую эффективность параллельного расчета поля интенсивности излучения методом дискретных направлений, причем применение технологии nVidia CUDA дает ускорение процедуры вычисления интеграла рассеяния в 82.3.
СПИСОК ЛИТЕРАТУРЫ
1. Черепащук А. М. SS 433: Новые результаты, новые проблемы // Земля и Вселенная. 1986. Т. 1. С. 21-29. [ A. M. Cherepashchuk, "SS 433: New results, new problems," (in Russian), Zemlya i Vselennaya, vol. 1, pp. 21-29, 1986. ]
2. A Summer of SS433: Forty Days of VLBA Imaging / Mioduszewski A. J. [et al.] // Bulletin of the American Astronomical Society. 2004. Vol. 36. P. 967. [ A. J. Mioduszewski, et al., "Forty Days of VLBA Imaging," Bulletin of the American Astronomical Society, vol. 36, pp. 967, 2004. ]
3. Галанин М. П., Лукин В. В., Чечеткин В. М. Ускорение джетов при различных вариантах моделирования источника вещества // Матем. моделирование. 2011. Т. 23, № 10. С. 65-81. [ M. P. Galanin, V. V. Lukin, and V. M. Chechetkin, "Jets acceleration in a variety of matter source modelling ways," Matem. Mod., vol. 23, no. 10, pp. 6581, 2011. ]
4. Savel'ev V. V., Toropin Yu. M., Chechetkin V. M.
A Possible Mechanism for the Formation of Molecular Flows // Astronomy Reports. 1996. Vol. 40. P. 494-508. [ M. P. Galanin, Yu. M. Toropin, and V. M. Chechetkin, "A Possible Mechanism for the Formation of Molecular Flows," Astronomy Reports, vol. 40, pp. 494-508, 1996. ]
5. Галанин М. П., Торопин Ю. М., Чечеткин В. М. Радиационное ускорение порций вещества в аккреционных воронках около астрофизических объектов // Астрономический журнал. 1999. Т. 76, № 2. С. 143-160. [ M. P. Galanin,
г
16
Yu. M. Toropin, and V. M. Chechetkin, "The Radiative Acceleration of Matter Bullets in Accretion Funnels near Astrophysical Objects", (in Russian), Astronomicheskii Journal, vol. 76, no. 2, pp. 143-160, 1999. ]
6. Четверушкин Б. Н. Математическое моделирование задач динамики излучающего газа. М.: Наука, 1985. 304 с. [ B. N. Chetverushkin, Mathematical modelling in radiative gas dynamics problems. Moscow: Nauka, 1985. ]
7. Pogorelov N. V., Semenov A. Yu. Solar wind interaction with the magnetized interstellar medium // Astron. Astrophys. 1997. Vol. 321. P. 330-337. [ N. V. Pogorelov and A. Yu. Semenov, "Solar wind interaction with the magnetized interstellar medium," Astron. Astrophys., vol. 321, pp. 330337, 1997. ]
8. 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. ]
9. Галанин М. П., Лукин В. В., Чечеткин В. М. Радиационное ускорение астрофизического канализированного струйного выброса // Вестник МГТУ им. Н. Э. Баумана. Естественные науки. Спец. вып. "Прикладная математика". 2011. С. 11-33. [ M. P. Galanin, V. V. Lukin, and V. M. Chechetkin, "Radiation Acceleration of Astrophysical Channelized Jet Outflows," (in Russian), Vestnik MGTU im. N. E. Baumana, Estestvennye nauki, SV-4, pp. 11-33, 2011. ]
10. Лукин В. В., Марчевский И. К. Учебно-экспериментальный вычислительный кластер. Часть 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. ]
ОБ АВТОРАХ
ГАЛАНИН Михаил Павлович, зав. отд. ИПМ, проф. каф. прикл. математики. Дипл. физик (МГУ им. М. В. Ломоносова, 1979). Д-р физ.-мат. наук по мат. моделированию (там же, 1995). Иссл. в обл. мат. моделир. и числ. методов.
ЛУКИН Владимир Владимирович, науч. сотр. ИПМ. Дипл. инж.-математик (МГТУ им. Н. Э. Баумана, 2008). Канд. физ.-мат. наук по мат. моделир. (ИПМ им. М. В. Келдыша, 2011). Иссл. в обл. мат. моделир. и числ. методов.
ЧЕЧЕТКИН Валерий Михайлович, гл. науч. сотр. ИПМ, проф. каф. теорет. ядерн. физики МИФИ. Д-р физ.-мат. наук, проф. Иссл. в обл. мат. моделир. астрофизических явлений и объектов.
METADATA
Title: Parallel software package for modelling of astrophysical jets using hybrid shared memory computation systems.
Authors: M. P. Galanin1, V. V. Lukin2, V. M. Chechetkin3.
Affiliation:
1,2,3 Keldysh Institute of Applied Mathematics RAS, Russia.
1,2 Bauman Moscow State Technical University (MGTU), Russia.
Email: 2vvlukin@gmail.com.
Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 18, no. 1 (62), pp. 230-240, 2014. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print).
Abstract: An RMHD model of the formation of astrophysical jet outflow from the vicinity of a compact object surrounded by a thin accretion disk and immersed into the nebula of galactic plasma is considered. The numerical algorithm is based on division by physical processes and includes a difference scheme for the MHD system of equations and discrete directions method for radiation transfer equation in the two-dimensional (2D) axisymmetrical case. Triangular unstructured staggered grids are used. Parallel numerical code using OpenMP and nVidia CUDA technologies for the shared memory systems is developed. The modeling results show the formation of the accelerating channel and acceleration of jet plasma up to 1/5 luminal speed and comply with the available observation data.
Key words: radiation magnetohydrodynamics; computations on GPU; astrophysical jets; shared memory systems.
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). CHECHETKIN, Valery Michailovich, Senior Researcher. KIAM RAS, Prof, Dept. of Theor. Nucl. Physics. Dr. of Phys.-Math. Sci.
