УДК 519.63, 524.3
AstroPhi: ПРОГРАММНЫЙ КОМПЛЕКС ДЛЯ МОДЕЛИРОВАНИЯ ДИНАМИКИ АСТРОФИЗИЧЕСКИХ ОБЪЕКТОВ НА ГИБРИДНЫХ СУПЕРЭВМ, ОСНАЩЕННЫХ УСКОРИТЕЛЯМИ Intel Xeon Phi1
И.М. Куликов, И.Г. Черных, Б.М. Глинский
В статье представлен новый сверхмасштабируемый программный комплекс AstroPhi для моделирования динамики астрофизических объектов на гибридных суперЭВМ, оснащенных ускорителями Intel Xeon Phi. Численный метод решения газодинамических уравнений основан на специально адаптированной для реализации на множестве ускорителей комбинации метода крупных частиц и метода Годунова. Для решения уравнения Пуассона используется быстрое преобразование Фурье. Программная реализация была отдельно протестирована на газодинамических задачах, на задаче решения уравнения Пуассона и на классических задачах гравитационной газовой динамики. Показано ускорение программного комплекса при использовании ускорителей Intel Xeon Phi, уточнено понятие масштабируемости при использовании ускорителей. Представлены результаты моделирования коллапса астрофизических объектов.
Ключевые слова: математическое моделирование, параллельные вычисления, ускорители Intel Xeon Phi, астрофизика.
Введение
Процессы коллапса астрофизических объектов в настоящее время активно исследуются теоретически в связи с появлением значительного числа наблюдательных данных. Явление коллапса имеет место как на начальной стадии звездной эволюции, так и на конечной стадии эволюции звезд (взрывы сверхновых с коллапсирующим ядром) [1]. Математическое моделирование играет более чем важную роль в теоретическом исследовании таких процессов. С каждым днем требования к астрофизическим моделям все возрастает и модели, которые несколько лет назад были актуальными, сейчас уже считаются устаревшими. Усложнение астрофизических моделей требует использования все больших вычислительных ресурсов, а, следовательно, модификации и создания новых вычислительных методов и параллельных алгоритмов для решения таких задач [2].
В последние два десятилетия из широкого диапазона газодинамических численных методов для решения нестационарных трехмерных астрофизических задач используются два основных подхода. Это лагранжев подход, в основном представленный SPH-методом [3, 4] (Smoothed Particle Hydrodynamics) и эйлеров подход с использование адаптивных сеток или AMR [5, 6] (Adaptive Mesh Refinement). В последние пять лет появился ряд программных пакетов с использованием комбинации лагранжева и эйлерова подходов. Основной проблемой SPH метода является поиск соседей и организация их гравитационного взаимодействия. Для эффективного решения этой задачи были разработаны ряд алгоритмов. Например, particle-particle/particle-mesh или P3M метод [7], адаптация P3M метода с использованием
1 Статья рекомендована к публикации программным комитетом Международной суперкомпьютерной кон-
ференции «Научный сервис в сети Интернет: все грани параллелизма — 2013».
иерархичности расчетной сетки AP3M [8], tree алгоритм [9], комбинация tree алгоритма и particle-mesh подхода Tree-PM метод [10]. Для решения уравнения Пуассона в сеточных методах используются в основном метод сопряженных градиентов (CGM), метод быстрого преобразования Фурье (FFT), метод последовательной верхней релаксации (SOR) и метод Федоренко [11] или многосеточный метод (MG).
Для численного решения газодинамических задач широкое применение получил метод Годунова [12], основным структурным элементом которого является задача о распаде произвольного разрыва (задача Римана) с параметрами газа в соседних ячейках разностной сетки. Как правило, параметры газа в соседних ячейках достаточно близки, что создает благоприятные условия для применения упрощенного алгоритма решения задачи о распаде разрыва. Различные алгоритмы получения приближенного решения задачи Римана дали большой класс методов [13, 14]. Основными методами являются методы типа Куранта-Изаксона-Риса [15] и Роу [16], которые строятся на основе использования различным образом линеаризованных гиперболических систем уравнений, Ошера [17], где решения задачи Римана строится только с использованием волн Римана.
Основоположным подходом к оценке скоростей волн является двухволновой метод Хартена-Лакса-Ван Лира (известный в литературе как HLL) [18], в котором учитываются левые и правые разрывы, без рассмотрения контактного разрыва. В схеме HLL использующей консервативные переменные, предложен простой, но эффективный способ выбора скоростей движения этих волн по максимальным наклонам характеристик в соседних ячейках разностной сетки. При этом веер волн разрежения заменяется скачком, но со скоростью распространения соответствующей максимальному наклону характеристик в этой волне разрежения. Поэтому этот выбор исключил проблему расчета «звуковой» точки при смене знака характеристик. Расчет ударных волн также проводится со скоростью превышающей точное значение. В этой связи схема HLL эффективна при расчете ударных волн и зон разрежения. Однако, в расчете энтропийных скачков методом установления принятое допущение приводит к неприемлемому «размазыванию» контактного разрыва. Существуют также модификации HLL, такие как HLLE [19], где особо учитываются крайние собственные числа линеаризованной задачи, и HLLC [20], где производится дополнительный учет центрального разрыва, движущегося со скоростью равной центральному собственному значению линеаризованной задачи Римана. Также на основе метода Годунова были сделаны реализации высокого порядка - монотонная противопотоковая схема второго порядка точности MUSCL [21] и TVD схемы [25], третьего порядка кусочно-параболический метод PPM [5]. Однако, что понимается под высоким порядком точности в случае разрывных решений, не очень понятно [26].
В ходе адаптации метода SPH к различным космологическим задачам возникло множество модификаций алгоритма. Для определения гидродинамических величин очень важен выбор радиуса сглаживания, который определяет количество влияющих на частицу соседей. Одним из свойств метода является то, что для получения верных решений необходимо сохранение числа соседей почти равным для всех частиц расчетной области, что не выполняется в случае коллапса. Если число соседей у всех элементов различается на несколько частиц, то результаты расчетов нельзя считать достоверными. Для преодоления такого недостатка метода вводится адаптивный радиус сглаживания, определяемый по числу соседей. Такое определение радиуса сглаживания приводит к ряду вычислительных сложностей, поэтому многие реализации метода допускают большие отклонения в числе со-
седей от частицы к частице. Таким образом, определение радиуса сглаживания допускает возможность выбора, а значит, оказывает влияние на решение. Авторы программ, использующие AMR, обычно задают величину наиболее подробного шага сетки по пространству, хотя сложно оценить степень необходимого сгущения сетки особенно в случае коллапса, поскольку в зависимости от степени адаптации сетки к решению в случае возникновения особенностей могут оставаться проблемы, характерные для сеток с линиями координат, располагающимися вдоль границ областей решения. Несомненно, разработка эйлеровых сеточных методов, не допускающих влияния сеточных линий на решение, позволит отказаться от методики AMR, а, значит, избежать всех вышеперечисленных недостатков этого подхода. Создание таких методов представляет собой более сложную задачу, чем построение адаптивных сеток, но, тем не менее, возможно. Как известно, уравнения газовой динамики инвариантны относительно некоторой группы точечных преобразований в пространстве независимых и зависимых переменных. Такая инвариантность является следствием инвариантности законов сохранений, из которых вытекают уравнения газовой динамики. Использование расчетной сетки неизбежно вносит неинвариантность в алгоритм расчета, что может оказывать влияние, например, на расчеты особенностей потока (ударных волн, контактных границ, слабых разрывов), которые движутся под различными углами к линиям сетки. Построению инвариантных относительно поворота разностных схем посвящены работы одного из авторов пакета [48, 49]
В рамках лагранжева подхода на основе SPH метода были разработаны пакеты Hydra [8], Gasoline [22], GrapeSPH [23], GADGET [24]. В рамках эйлерова подхода (в том числе и с использованием AMR) были разработаны пакеты NIRVANA [27], FLASH [28], ZEUS-MP [29], ENZO [6], RAMSES [30], ART [31], Athena [32], Pencil Code [33], Heracles [39], Orion [40], Pluto [41], CASTRO [42]. Эйлеров подход с использованием AMR был впервые использован на гибридных суперкомпьютерах, оснащенных графическими ускорителями, в пакете GAMER [34]. Пакет BETHE-Hydro [35], AREPO [36], CHIMERA [37] и авторский пакет PEGAS [38] основаны на комбинации лагранжева и эйлерова подходов.
Большое количество программных реализаций и численных методов говорит об актуальности исследований в области разработки новых методов и их программных реализаций для решения задач астрофизики. Кроме того, несмотря на развитие астрофизических пакетов в сторону петафлопсных вычислений, таких как PetaGADGET [43], Enzo-P [44], PetaART [45] нужно отметить фундаментальные ограничения в масштабируемости AMR и SPH подходов, которые используются в основном числе программных пакетов для решения задач астрофизики [46, 47].
Целью данной статьи является описание модификации численного метода и особенностей реализации пакета AstroPhi на гибридных суперкомпьютерах, оснащенных ускорителями Intel Xeon Phi, исследование его масштабируемости и исследование задач коллапса астрофизических объектов. Статья организована следующим образом. В разделе 1 описана численная схема решения уравнений гравитационной газовой динамики. Раздел 2 содержит описание параллельной реализации численной схемы. В разделе 3 приведена верификация программной реализации. Раздел 4 посвящен вычислительным экспериментам по иодели-рованию процесса коллапса астрофизических объектов. В заключении приведены выводы по созданному программному пакету и обозначены направления дальнейшего его развития.
1. Описание численной схемы
Будем рассматривать трехмерную модель динамики самогравитирующего газа в декартовых координатах, включающих в себя расширенную систему уравнений газовой динамики в дивергентной форме, замкнутую уравнением состояния для идеального газа. Система уравнений газовой динамики дополнена уравнением Пуассона для гравитационного потенциала и вкладом в потенциал от центрального тела. Уравнения записаны в безразмерном виде
где р — давление, р — плотность, V — вектор скорости, рЕ — плотность полной энергии, Ф — собственный гравитационный потенциал, Фо — вклад в гравитационный потенциал от центрального тела, е — внутренняя энергия, 7 — показатель адиабаты. В качестве основных характерных параметров выбраны радиус солнца Ь = Я0, масса солнца Мо = М&, гравитационная постоянная О = 6, 67 ■ 10-11 Н-м2/кг.
1.1. Метод решения уравнений газовой динамики
Введем в трехмерной области решения равномерную прямоугольную сетку с ячейками
Х% = ^hx, % = 1, ■■, ^шах) ук = кНу, к = 1, .., Кшах) %1 = , I = 1, .., ЬШах) где hx, Ну, hz — шаги
сетки, 1шах, Кшах, Ьшах - количество узлов сетки по направлениям х, у, г: Нх = Хшах/1шах, Ну = ушах/Кшах, hz = гшах/Ьшах. За основу метода решения системы уравнений газовой динамики выбран метод крупных частиц [38], уже хорошо зарекомендовавший себя в ходе решения астрофизических задач [2].
Исходная система газодинамических уравнений решается в два этапа. Система уравнений на первом, эйлеровом, этапе описывает процесс изменения параметров газа в произвольной области течения за счет работы сил давления, а также за счет разности потенциалов и охлаждения. Для исключения влияния направлений координатных линий использованы элементы операторного подхода [48]. В его основе определение плотности, давления, потенциала и импульса в ячейках, в узлах ячеек размещаются только вектор скорости. Для дискретных аналогов компонент скорости, определеных в узлах сетки, применяется функция осреднения в ячейку. Значения давления и скорости на всех границах ячеек Р и V суть точное решение линеаризованной системы уравнений эйлерова этапа по каждому из направлений осей координат без учета вклада потенциала и охлаждения:
йт(рЕУ) = —йіу(ру) — (рдтай(Ф + Ф0),У),
дтай(р) — рдтай(Ф + Фо),
ду 1 др
ді рдп
Эта система на каждой границе ячеек уже является линейной гиперболической системой и имеет аналитическое решение:
V _ УЬ + У К + РЬ - РЕ
Рь + Ре
РЬРЕ Ь - 1)(РЬ + РЕ)’
р _ РЬ + РЕ + УЬ - УЕ рь РЕ (1 - 1)(РЬ + РЕ)
2 2 у рь + РЕ ’
где /ь, /е — значения соответствующих функций (р, У,р) справа и слева от границы ячеек. Эти значения и используются в схеме эйлерового этапа.
Система уравнений на втором, лагранжевом, этапе, содержит дивергентные слагаемые вида
д. + (1у(/У) _ 0
и отвечает за процесс адвективного переноса всех газодинамических величин /. В исходном варианте численного метода использовался подход, связанный с вычислением вклада в соседние ячейки со схемной скоростью [49]. Однако, такой подход совершенно не годится в случае использования ускорителей. Для этого рассмотрим решение следующей одномерной постановки предыдущего уравнения:
ги+1 Ги
■) ікі г ікі
Т
ри+1/2 ^п+1/2
+ і+1/2,кі і-1/2,кІ о
к
где величина определяется по правилу полной деформации ячейки (см. рис. 1):
К
1+1/2 _ ^ Уі+1/2,к±1,1±11,
+
ікі
і+1/2,кі
Іікі,Уі+1/2,к±1,і±1 > 0 Іі+1,кі,Уі+1/2,к±1,і±1 < 0-
Рис. 1. Поток соответствующий газодинамической величины
На каждом временном шаге производится корректировка баланса энергий [50]. С этой целью осуществляется перенормировка схемных скоростей переноса массы, импульса и двух видов энергий на лагранжевом этапе метода таким образом, что происходит корректировка длины вектора скорости при неизменном направлении. Для этого все компоненты скорости в каждой ячейке области умножаются на множитель ^,/,1:
КікІ —
\
х,ікі
+ у2у,ікІ + У
2
г'ікі
4
2
Такая модификация метода обеспечивает справедливость детального баланса энергий. Заметим, что разностная схема не становится полностью консервативной, поскольку коррекция скорости вносит погрешность в закон сохранения импульса. Модификация схемы позволяет производить расчет одного шага по времени независимо для всех ячеек расчетной области. Таким образом, будет полностью отсутствовать взаимодействие между ядрами ускорителя, что должно привести к большому ускорению при использовании Intel Xeon Phi.
1.2. Метод решения уравнения Пуассона
После реализации газодинамической системы уравнений решается уравнение Пуассона для гравитационного потенциала. Для его решения используется 27-точечный шаблон. Потенциал и плотность представляется в виде суперпозиции по собственным функциям оператора Лапласа. Получим следующую схему решения уравнения Пуассона в пространстве гармоник:
ф ^ pjmn
3mn= 6(1 - (1 - 3 sin2(^))(1 - 3 sin2(K))(1 - | sin2(^))) •
Таким образом схема решения уравнения Пуассона примет следующий вид:
1. Преобразование в пространство гармоник выражения.
2. Решение в пространстве гармоник уравнения.
3. Обратное преобразование из пространства гармоник функции потенциала.
Для перехода в пространство гармоник и обратно, которое состоит в нахождении коэффициентов перехода, воспользуемся быстрым преобразованием Фурье.
Краевые условия уравнения Пуассона определяют решение задачи, поэтому их постановка является достаточно важной проблемой. Известно, что в бесконечном удалении от объекта гравитационный потенциал может считаться нулевым. Краевые условия приходится ставить на конечном расстоянии от газового объекта. Для решения данной задачи был предложен следующий вариант: считать, что масса тела сосредоточена в центре рассматриваемой области и модуль потенциала обратно пропорционален расстоянию от рассматриваемой границы до центра области.
2. Описание параллельной реализации
Трехмерность модели и нестационарность задачи выдвигают строгие требования к экономичности используемых методов решения. В последнее время бурное развитие вычислительной техники позволило производить ресурсоемкие расчеты и получать физически оправданные результаты для трехмерных программ. Использование суперкомпьютеров позволяет использовать большие объемы данных, на порядки повышать производительность вычислений, а как следствие, и точность. Основные вычислительные затраты приходятся на решение гидродинамических уравнений, решение которых занимает 90 процентов времени (рис. 2).
2.1. Параллельная реализация гидродинамических уравнений
В основе параллельной реализации решения гидродинамических уравнений лежит многоуровневая одномерная декомпозиция расчетной области. По одной координате внешнее
Рис. 2. Процентное соотношение вычислительных затрат на решение каждого из этапов
одномерное разрезание происходит средствами технологии МР1, внутри каждой подобласти разрезание происходит средствами ОрепМР, адаптированного для М1С-архитектур (рис. 3).
Рис. 3. Схема декомпозиции расчетной области для решения гидродинамических уравнений
Это связано с топологией и архитектурой гибридного суперЭВМ МВС-10П межведомственного суперкомпьютерного центра РАН, который был использован для вычислительных экспериментов. Модификация численного метода решения гидродинамических уравнений позволяет на каждом этапе численного метода независимо вычислять значения потоков через каждую ячейку. Декомпозиция области на каждом этапе осуществляется с перекрытием одного слоя граничных точек соседних областей.
2.2. Параллельная реализация уравнения Пуассона
Трехмерное параллельное быстрое преобразование Фурье выполняется с помощью процедуры из свободно распространяемой библиотеки FFTW. Способ распределения массивов также задается библиотекой. Перекрытие расчетных областей не требуется. В силу малых вычислительных затрат на решение уравнения Пуассона по сравнению с решением гидродинамических уравнений ускорители не использовались (на рис. 4 показано ускорение решения уравнения Пуассона для расчетной сетки 10243 в зависимости от числа процессов). Однако, в дальнейшем такая реализация, основанная на архитектуре библиотеки FFTW и адаптированная для решения задачи, обязательно будет сделана.
2.3. Эффективность параллельной реализации
Для определения эффективности гибридной реализации мы вводим следующие понятия масштабируемости. SingleMIC performance (сильная масштабируемость в рамках одного ускорителя Intel Xeon Phi) - уменьшение времени счета одного шага одной и той же
Рис. 4. Ускорение решения уравнения Пуассона
задачи при использовании большего числа ядер ускорителя. MultiMIC performance (слабая масштабируемость при использовании многих ускорителей Intel Xeon Phi) - сохранение времени счета одного шага одного и того же объема задачи при одновременном увели-ченпп количества ускорителей. FFTW performance (сильная масштабируемость при использовании библиотеки FFTW) - уменьшение времени счета одного шага одной п той же задачи при использовании большего числа процессоров или ядер. Результаты эффективности программной реализации приведена на рис. Б (на рисунке изображены эффективность параллельной реализации газодинамического этапа в зависимости от использованных ускорителей (а) п ускорение на одном Intel Xeon Phi для расчетной сеткп 1283 (б)).
1,00-
.о
I- 0,9В -\
О
о
СО 0,96 -
(D 0,94
-80 0,92-1
0,90
ч
6 12 18 24
Число ускорителей
а)
30
30-|
24-
ш
т 18-
ф
LL
О
VS-
О
>
6-
0-
/
Г
/
10
20 30 40
Число ядер
б)
50 ео
Рис. 5. Эффективность п ускорение решения газодинамических уравнений
Таким образом, для самого тяжелого - газоднамического этапа численного метода -было получено двадцатисемикратное ускорение при полном использовании ускорителя Intel Xeon Phi.
3. Верификация программной реализации
Так как в астрофизике математическое моделирование зачастую выступает единственной возможностью подтвердить или опровергнуть новые теории, то исследователи особен-
64 Вестник ЮУрГУ. Серия «Вычислительная математика и информатика»
но нуждаются в применении надежных и заслуживающих доверия программ. Прежде чем представлять новые результаты моделирования, необходимо провести разнообразные тестовые расчеты для обоснования и верификации используемой программы. Верификация и обоснование - основные этапы развития для любой технологии, будь это пакет программ для математического моделирования или инструментарий для наблюдений. Для вычислительной технологии целью такого этапа тестирования является оценка правомерности и точности моделирования.
3.1. Тесты Годунова
В области вычислительной гидродинамики проделана большая работа по обоснованию и верификации. В процессе создания комплекса программ проводилась верификация численного алгоритма на тестах с решениями из специализированного банка данных. Одна из проблем - моделирование ударной волны. Известно, что различные методы по-разному моделируют эту область, а именно либо с осцилляциями, либо с диссипацией. В случае гравитационной газовой динамики осцилляции являются менее желательными, так как любая, в общем случае, случайная волна плотности будет стягивать на себя соседний газ, что приведет к образованию нефизичных флуктуаций газодинамических параметров. Также существует проблема моделирования существенной области разрежения. Известно, что многие методы дают нефизичный рост внутренней энергии в этом месте. Кроме этого стандартной проверкой робастности метода служит огромный начальный перепад давления (пять десятичных порядков), который должен выявить способность метода устойчиво моделировать сильные возмущения с возникновением быстро распространяющихся ударных волн. На всех этих тестах (см. табл. 1) была верифицирована данная программная реализация.
Таблица
Начальная конфигурация ударной трубы
Параметр Тест 1 Тест 2 Тест 3
рь 2 1 1
Уь 0 -2 0
Рь 2 0,4 1000
ря 1 1 1
Уя 0 2 0
Ря 1 0,4 0,01
Хо 0,5 0,5 0,5
і 0,2 0,15 0,012
Целью первого теста является определение правильности описания контактного разрыва. Большинство методов решения газодинамических уравнений дают либо осцилляцию, либо диффузию («размазывание» ударных волн). Авторский метод дает размазывание решения в области контактного разрыва, которое уменьшается с дроблением сетки. В ходе второго теста, газ с одинаковыми термодинамическими параметрами разлетается в разные стороны, образуя в центре существенную область разрежения. Тест выявляет способность физически правдоподобно моделировать такую ситуацию. Из литературы известно, что
многие методы дают ошибочный (нефизический) рост температуры в области сильного разрежения и как следствие, получаемое решение искажается. Авторский метод успешно моделирует область разрежения. Основная задача третьего теста - проверка устойчивости численного метода. Огромный перепад давления (5 десятичных порядков) должен выявить способность метода устойчиво моделировать сильные возмущения с возникновением быстро распространяющихся ударных волн. Результаты показывают, что имеют место малые осцилляции решения в области контактного разрыва. Так называемая волна-предшественник (ступенька на графике внутренней энергии на правом фронте ударной волны) отражена корректно, без размазывания, что говорит в пользу метода (см. рис. 6).
2,0
1,8-
1,6-
1,4-
1,2-
1,0-
0,0 0,2 0,4 0,6 0,8 1,0
а)
б)
0,0 0,2 0,4 0,6 0,8 1,0
в)
2,0
Фгсоеееоаавеаеап
0,0 0,2 0,4 0,6 0,8 1,0
д)
е)
0,3
0,0
25-
1,8- ■
1,2- ЛУ 20"
/ о°° °° ■
/о 0 0,6- 15-
0
° 0 0,0- ■
/> ( 10-
-0,6- ■
7 -1,2- 5
о° о .
° ° -18-
,п^п0о° 1
0,0 0,2 0,4 0,6 0,8 1,0
ж)
0,0 0,2 0,4 0,6
з)
0,6 0,8 1,0 Координата
и)
Рис. 6. Численное решение тестов Годунова
На рисунке 6 показаны распределения плотности, скорости и давления в результате моделирования первого теста (а, г, ж соответственно), второго теста (б, д, з соответственно), третьего теста (в, е, и соответственно), сплошной линией обозначено точное решение, точками обозначен результат расчета.
3.2. Тест Аксенова
Рассмотрим систему уравнений одномерной газовой динамики в размерном виде:
ди ди 1 др
дЬ + и дх рдх0
др + дри = 0 дЬ + дх 0
p
Ро \Ро)
где р - давление, р - плотность, и - скорость, 7 - показатель адиабаты.
В качестве характерных величин выберем I - характерная длина, ро - характерная
подхода, описанного в [51], можно взять в качестве начальных данных r = 1 + 0, 5 cos x, z = 0. Тогда периодическое решение на интервале [0; 2п] записывается в виде:
Легко проверить, что такое решение с учетом обезразмеривания удовлетворяет исходной системе уравнений. Для сравнения численного результата, полученного авторским методом, с аналитическим решением, выберем момент времени, когда хотя бы одна из функций имеет простой (явный) вид. Выберем в качестве такого момента £ = п/2. Тогда г(х) = 1, а уравнение для г имеет простой вид:
Результаты вычислительного эксперимента представлены на рис. 7 (на рисунке изображены распределения плотности (а) и скорости (б) в момент времени £ = п/2, сплошной линией обозначено точное решение, точками обозначен результат расчета).
Из рисунков видно, что решение скорости ввиду его гладкости достаточно хорошо приближается численным решением, а график плотности имеет скачок в центре. Данный скачок имеет ту же природу, что и скачок температуры в третьем тесте Годунова, когда газ разлетается в разные стороны. Фактически в этом тесте имеет место энтропийный след, который образуется в результате «стекания» газа в эту область с нулевой скоростью. Полученную особенность имеет смысл исследовать для сравнения численных методов на таком тесте и возникающих особенностей при воспроизведении данной области. В любом случае стоит отметить, что данная особенность сосредоточена на конечном числе точек расчетной области и уменьшается при дроблении сетки.
плотность, р0 - характерное давление. В этом случае характерная скорость и0 = у/7р0/ро,
и характерное время £0 = I/у/7р0/р0. Выбрав в качестве размерных величин I = 1, р0 = 1, р0 = 1, 7 = 3 и введя обозначения Л = 1/(7 — 1), г = р1/2Л, г = и/2 Л, с использованием
r = І + О, Б cos(x — zt) cos(rt),
z = О, Б sin(x — zt) sin(rt).
z = 0,5sin(x — zt).
0,6-1
[
-0,6
0
2 3 4 5 6
Координата
0
2 3 4 5
Координата
6
б)
Рис. Т. Численное решение теста Аксенова
3.3. Неустойчивость Кельвина-Гельмгольца и Релея—Тейлора
В основе физической постановки задачи лежит гравитационная неустойчивость, что приводит к математической некорректности по Адамару. Численный метод не может подавлять физическую неустойчивость. Для проверки корректного воспроизведения неустойчивых течений была сделана верификация численного метода на задачах о развитиии неустойчивости Релея-Тейлора и Кельвина-Гельмгольца. В случае моделирования неустойчивости Релея-Тейлора проверяется возможность воспроизведения гравитационного терма. Неустойчивость Кельвина-Гельмгольца позволяет убедится в возможности метода моделировать нелинейную гидродинамическую турбулентность.
Начальные условия для моделирования неустойчивости Релея-Тейлора: [—0, 5; 0, 5]2 -область моделирования, 7 = 1, 4 - показатель адиабаты,
ро(х) = ( 1 Г - 0 | 2, г > 0
р = 2, 5—рду - равновесное давление, д - ускорение свободного падения, %0(х, у) = А(у)[1 + 008(2пх)][1 +008(2пу)], где
Л10-2, М< 01,
\ 0, у> 0,01
Начальные условия для моделирования неустойчивости Кельвина-Гельмгольца:
[—0, 5; 0, 5]2 - область моделирования, 7 = 1, 4 - показатель адиабаты,
Ро(х) = ( 1 Г < 0, о 2, г > 0
„ = / 0, 5, |у| < 0,25,
х | —0, 5, |у| > 0, 25
р = 2, 5 - начальное давление, %0(х,у) = А(у)[1 + еов(8пх)][1 + ео8(8пу)], где
А(у) = / 10-2, ||у|— 0, 251 — 0, 01,
\ 0, ||у|— 0, 25| > 0, 01
Результаты моделирования неустойчивости Кельвина-Гельмгольца и Релея-Тейлора представлены на рис. 8.
4. Моделирование процесса коллапса астрофизических объектов
Процессы коллапса астрофизических объектов в настоящее время активно исследуются теоретически в связи с появлением значительного числа наблюдательных данных. Явление коллапса имеет место как на начальной стадии звездной эволюции, так и на конечной стадии эволюции звезд (взрывы сверхновых с коллапсирующим ядром).
4.1. Сжатие невращающегося облака
Приведем результаты вычислительного эксперимента моделирования коллапса. Сначала необходимо оценить точность моделирования коллапса. Основной критерий правильно-
68 Вестник ЮУрГУ. Серия «Вычислительная математика и информатика»
-0,4 -0,2 0,0 0,2 0,4 -0,4 -0,2 0,0 0,2 0,4
а) б)
Рис. 8. Моделирование неустойчивости Кельвина-Гельмгольца (а) и Релея-Тейлора (б)
сти - поведение полной энергии системы. Затем необходимо провести сравнение профилей плотности, полученных с помощью реализаций методов FLIC (авторская реализация метода) и SPH. Профиль начального распределения плотности р (r) ~ 1/r. Начальное распределение давления p = 0,1 • р7, показатель адиабаты y = 5/3.
Целью исследования точности моделирования коллапса является поведение полной энергии при дроблении сетки. Источником ошибки в законе сохранения полной энергии является конечная стадия коллапса, когда плотность и другие газодинамические величины увеличиваются в 10 или 100 раз. Подавляющая масса газа находится в шаре с радиусом Rcollaps = 0, 1 • Ro (10 % от начального радиуса газового шара). Поэтому для моделирования этого процесса нам нужно иметь на радиусе Rcoiiaps достаточное число ячеек для удовлетворительного моделирования процесса. Проведем сравнение поведения закона сохранения полной энергии при дроблении сетки (см. рис. 9).
Рис. 9. Относительная погрешность полной энергии
Из рисунка видно, что при дроблении сетки относительная погрешность полной энергии уменьшается, а, значит, дробление сетки приведет к моделированию коллапса с наперед заданной точностью. На сетке 512 х 512 х 512 погрешность на уровне 5 %, что является удовлетворительным для качественного сравнения решения полученного FLIC методом с решением, полученным SPH методом.
AstroPhi: программный комплекс для моделирования динамики астрофизических...
4.2. Сжатие быстро вращающегося облака
В рамках исследования возможности моделирования коллапса вращающихся прото-звездных облаков будем моделировать газовое облако, ограниченное сферой радиуса Е0 =
3, 81 ■ 1014 м, с массой Мд = 3, 457 ■ 1030 кг, с равномерным распределением плотности р = 1, 492 ■ 10-14 кг/м3 и давления р = 0,1548 ■ 10-10 Н/м2, вращающийся с угловой скоростью ш = 2, 008 ■ 10-12 рад/с. Показатель адиабаты соответствует водороду 7 = 5/3. Масса центрального тела М0 = 1, 998 ■ 1030 кг. В качестве размерных величин выберем следующие значения: Ь0 = 3, 81 ■ 1014 м, р0 = 1, 492 ■ 10-14 кг/м3, р = 0,1548 ■ 10-7 Н/м2, -и0 = 1010 м/с, £ = 3, 7 ■ 1011 с, ш0 = 0, 27 ■ 10-11 рад/с. Тогда в безразмерных величинах задача ставится следующим образом: р = 1, 0 - плотность газового облака, р = 10-3 - давление в газовом облаке, ш = 0, 744 - угловая скорость вращения, шС0 = 2, 42 - масса центрального тела, 7 = 5/3 - показатель адиабаты, [0; 6, 4]3 - расчетная область.
В ряде работ, посвященных коллапсированию протозвездных облаков, проводился поиск ответа, каким будет распределение плотности в экваториальной плоскости облака после коллапсирования. Четкого ответа так и не было получено. Основным результатом в данной задаче является поведение энергий (см. рис. 10) и торобразный профиль плотности в экваториальной плоскости.
0,0 0,5 1,0 1,5 2,0
Время
Рис. 10. Поведение различных видов энергии в задаче сжатия быстро вращающегося облака
Стоит отметить, что поведение энергий качественно, а до момента коллапса - количественно, - совпадают с результатами других авторов [1].
4.3. Сжатие вращающегося молекулярного облака
В рамках исследования возможности моделирования коллапса вращающихся молекулярных облаков будем моделировать газовое облако, ограниченное сферой радиуса Е0 = 100 парсек, с массой Мд = 107М0, с распределением плотности р (г) ~ 1/г и температуры Т ~ 2000 К, вращающийся с угловой скоростью ш = 21 км/с. Показатель адиабаты соответствует водороду 7 = 5/3. Скорость звука с ~ 3, 8 км/с. В качестве размерных величин выберем следующие значения: £0 = 100 парсек, р0 = 1, 2 ■ 10-18 кг/м3, -и0 = 21 км/с. Тогда в безразмерных величинах задача ставится следующим образом: р = 1, 0 - плотность газового облака в центре, р = 2 х 10-2 - давление в газовом облаке в центре, ш = 1 - угловая скорость вращения, 7 = 5/3 - показатель адиабаты, [0;6, 4]3 - расчетная область.
0,0 0,5 1,0 1,5 2,0 2,5
Время
Рис. 11. Поведение различных видов энергии в задаче сжатия молекулярного облака
В рамках данного исследования количественно поведение энергий (см. рис. 11) совпало с результатом других авторов [53].
Заключение
Статья посвящена новому сверхмасштабируемому программному комплексу AstroPhi для моделирования динамики астрофизических объектов на гибридных суперЭВМ, оснащенных ускорителями Intel Xeon Phi. В статье описан новый численный метод решения газодинамических уравнений, основанный на специально адаптированной для реализации на множестве ускорителей комбинации метода крупных частиц и метода Годунова. Для решения уравнения Пуассона используется быстрое преобразование Фурье. Программная реализация была отдельно протестирована на газодинамических задачах, на задаче решения уравнения Пуассона и на классических задачах гравитационной газовой динамики. Показано ускорение программного комплекса при использовании ускорителей Intel Xeon Phi, уточнено понятие масштабируемости при использовании ускорителей. Представлены результаты моделирования коллапса астрофизических объектов. В текущей реализации AstroPhi реализована только модель гравитационной газовой динамики, в дальнейшем планируется добавить магнитную составляющую и бесстолкновительную компоненту, основанную на первых моментах уравнения Больцмана.
Исследование выполнено при частичной финансовой поддержке грантов РФФИ № 1201-31352 мол-а, 13-07-00589-а, 12-07-00065-а, 13-01-00231-а, 11-05-00937; МИП № 39 СО РАН, МИП № 130 СО РАН; гранта Президента РФ МК-4183.2013.9; Программы Президиума РАН Проект № 15.9; гранта Мэрии города Новосибирска.
Литература
1. Ardeljan, N.V. An Implicit Lagrangian Code for the Treatment of Nonstationary Problems in Rotating Astrophysical Bodies / N.V. Ardeljan, G.S. Bisnovatyi-Kogan, S.G. Moiseenko // Astronomy & Astrophysics. — 1996. — Vol. 115. — P. 573-594.
2. Tutukov, A. Gas Dynamics of a Central Collision of Two Galaxies: Merger, Disruption, Passage, and the Formation of a New Galaxy / A. Tutukov, G. Lazareva, I. Kulikov //
Astronomy Reports. — 2011. — Vol. 55, I. 9. — P. 770-783.
3. Gingold, R.A. Smoothed Particle Hydrodynamics: Theory and Application to Non-spherical Stars / R.A. Gingold, J.J. Monaghan // Monthly Notices of the Royal Astronomical Society.
— 1977. — Vol. 181. — P. 375-389.
4. Luci, L.B. A Numerical Approach to the Testing of the Fission Hypothesis / L.B. Luci // The Astrophysical Journal. — 1977. — Vol. 82, № 12. — P. 1013-1024.
5. Collela, P. The Piecewise Parabolic Method (PPM) for Gas-dynamical Simulations / P. Collela, P.R. Woodward // Journal of Computational Physics. — 1984. — Vol. 54. — P. 174-201.
6. Norman, M. The Impact of AMR in Numerical Astrophysics and Cosmology / M. Norman // Lecture Notes in Computational Science and Engineering. — 2005. — Vol. 41. — P. 413-430.
7. Hockney, R.W. Computer Simulation Using Particles / R.W. Hockney, J.W. Eastwood — New York: McGraw-Hill, 1981. — 540 p.
8. Couchman, H. Hydra Code Release / H. Couchman, F. Pearce, P. Thomas. URL: http://arxiv.org/abs/astro-ph/9603116 (дата обращения: 25.07.2013).
9. Barnes, J. A Hierarchical O(N log N) Force-calculation Algorithm / J. Barnes, P. Hut // Nature. — 1986. — Vol. 324. — P. 446-449.
10. Dubinski, J. GOTPM: A Parallel Hybrid Particle-Mesh Treecode / J. Dubinski, J. Kim, C. Park, R. Humble // New Astronomy. — 2004. — Vol. 9. — P. 111-126.
11. Fedorenko, R. A Relaxation Method for Solving Elliptic Difference Equations / R. Fedorenko // U.S.S.R. Computational Mathematics and Mathematical Physics. — 1961. — Vol. 1. — P. 1092-1096.
12. Godunov, S.K. A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations / S.K. Godunov // Matematicheskii Sbornik. — 1959. — Vol. 47.
— P. 271-306.
13. Kulikovskii, A.G. Mathematical Aspects of Numerical Solution of Hyperbolic Systems /
A. Kulikovskii, N. Pogorelov, A. Semenov. — Moscow:Fizmatlit, 2001. — 608 p.
14. Toro, E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics / E. Toro. — Heidelberg:Springer-Verlag, 1999. — 686 p.
15. Courant, R. On the Solution of Nonlinear Hyperbolic Differential Equations by Finite Difference / R. Courant, E. Isaacson, M. Rees // Communications on Pure and Applied Mathematics. — 1952. — Vol. 5, № 3. — P. 243-255.
16. Roe, P. Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes / P. Roe // Journal of Computational Physics. — 1997. — Vol. 135, I. 2. — P. 250-258.
17. Engquist, B. One-sided Difference Approximations for Nonlinear Conservation Laws /
B. Engquist, S.J. Osher // Mathematics of Computation. — 1981. — Vol. 36, № 154. — P. 321-351.
18. Harten, A. On Upstream Differencing and Godunov-type Schemes for Hyperbolic Conservation Laws / A. Harten, P.D. Lax, B. Van Leer // Society for Industrial and Applied Mathematics Review. — 1983. — Vol. 25, № 1. — P. 35-61.
19. Einfeld, B. On Godunov-type Methods for Gas Dynamics / B. Einfeld // Society for Industrial and Applied Mathematics Journal on Numerical Analysis. — 1988. — Vol. 25, № 2. — P. 294-318.
20. Batten, P. On the Choice of Savespeeds for the HLLC Riemann Solver / P. Batten, N. Clarke,
C. Lambert, D.M. Causon // Society for Industrial and Applied Mathematics Journal on Computing. — 1997. — Vol. 18, № 6. — P. 1553-1570.
21. Van Leer, B. Towards the Ultimate Conservative Difference Scheme. V - A Second-order Sequel to Godunov’s Method / B. Van Leer // Journal of Computational Physics. — 1979.
— Vol. 32. — P. 101-136.
22. Wadsley, J.W. Gasoline: a Flexible, Parallel Implementation of TreeSPH / J.W. Wadsley, J. Stadel, T. Quinn // New Astronomy. — 2004. — Vol. 9, I. 2. — P. 137-158.
23. Matthias, S. GRAPESPH: Cosmological Smoothed Particle Hydrodynamics Simulations with the Special-purpose Hardware GRAPE / S. Matthias // Monthly Notices of the Royal Astronomical Society. — 1996. — Vol. 278, I. 4. — P. 1005-1017.
24. Springel, V. The Cosmological Simulation Code GADGET-2 / V. Springel // Monthly Notices of the Royal Astronomical Society. — 2005. — Vol. 364, I. 4. — P. 1105-1134.
25. Jin, S. The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions / S. Jin, Z. Xin // Communications on Pure and Applied Mathematics. — 1995. — Vol. 48. — P. 235-276.
26. Godunov, S.K. Experimental Analysis of Convergence of the Numerical Solution to a Generalized Solution in Fluid Dynamics / S.K. Godunov, Yu.D. Manuzina, M.A. Nazareva // Computational Mathematics and Mathematical Physics. — 2011. — Vol. 51. — P. 88-95.
27. Ziegler, U. Self-gravitational Adaptive Mesh Magnetohydrodynamics with the NIRVANA Code / U. Ziegler // Astronomy & Astrophysics. — 2005. — Vol. 435. — P. 385-395.
28. Mignone, A. The Piecewise Parabolic Method for Multidimensional Relativistic Fluid Dynamics / A. Mignone, T. Plewa, G. Bodo // The Astrophysical Journal. — 2005. — Vol. 160. — P. 199-219.
29. Hayes, J. Simulating Radiating and Magnetized Flows in Multiple Dimensions with ZEUS-MP / J. Hayes, et al. // The Astrophysical Journal Supplement Series. — 2006. — Vol. 165.
— P. 188-228.
30. Teyssier, R. Cosmological Hydrodynamics with Adaptive Mesh Refinement. A New High Resolution Code Called RAMSES / R. Teyssier // Astronomy & Astrophysics. — 2002. — Vol. 385. — P. 337-364.
31. Kravtsov, A. Constrained Simulations of the Real Universe. II. Observational Signatures of Intergalactic Gas in the Local Supercluster Region / A. Kravtsov, A. Klypin, Y. Hoffman // The Astrophysical Journal. — 2002. — Vol. 571. — P. 563-575.
32. Stone, J. Athena: A New Code for Astrophysical MHD / J. Stone, et al. // The Astrophysical Journal Supplement Series. — 2008. — Vol. 178. — P. 137-177.
33. Brandenburg, A. Hydromagnetic Turbulence in Computer Simulations / A. Brandenburg, W. Dobler // Computer Physics Communications. — 2002. — Vol. 147. — P. 471-475.
34. Schive, H. GAMER: a GPU-accelerated Adaptive-Mesh-Refinement Code for Astrophysics j H. Schive, Y. Tsai, T. Chiueh jj The Astrophysical Journal. — 20І0. — Vol. 18б. — P. 457-484.
35. Murphy, J. BETHE-Hydro: An Arbitrary Lagrangian-Eulerian Multidimensional
Hydrodynamics Code for Astrophysical Simulations j J. Murphy, A. Burrows jj The Astrophysical Journal Supplement Series. — 2008. — Vol. 179. — P. 209-241.
36. Springel, V. E Pur Si Muove: Galilean-invariant Cosmological Hydrodynamical Simulations on a Moving Mesh j V. Springel jj Monthly Notices of the Royal Astronomical Society. —
2010. — Vol. 401. — P. 791-851.
37. Bruenn, S. 2D and 3D Core-collapse Supernovae Simulation Results Obtained with the CHIMERA Code j S. Bruenn, et al. jj Journal of Physics. — 2009. — Vol. 180. — P. 1-5.
38. Vshivkov, V. Hydrodynamical Code for Numerical Simulation of the Gas Components of Colliding Galaxies j V. Vshivkov, G. Lazareva, A. Snytnikov, I. Kulikov, A. Tutukov jj The Astrophysical Journal Supplement Series. — 2011. — Vol. 194, I. 47. — P. 1-12.
39. Gonzalez, M. HERACLES: a Three-dimensional Radiation Hydrodynamics Code j M. Gonzalez, E. Audit, P. Huynh jj Astronomy к Astrophysics. — 2007. — Vol. 4б4. — P. 429-435.
40. Krumholz, M.R. Radiation-hydrodynamic Simulations of the Formation of Orion-like Star Clusters. I. Implications for the Origin of the Initial Mass Function j M.R. Krumholz, R.I. Klein, C.F. McKee, J. Bolstad jj The Astrophysical Journal. — 2007. — Vol. бб7, I. 74. — P. 1-1б.
41. Mignone A. PLUTO: A Numerical Code for Computational Astrophysics j A. Mignone, et al. jj The Astrophysical Journal Supplement Series. — 2007. — Vol. 170. — P. 228-242.
42. Almgren A. CASTRO: A New Compressible Astrophysical Solver. I. Hydrodynamics and Self-gravity j A. Almgren, et al. jj The Astrophysical Journal. — 2010. — Vol. 715. — P. 1221-1238.
43. Feng Y. Terapixel Imaging of Cosmological Simulations j Y. Feng, et al. jj The Astrophysical Journal Supplement Series. — 2011. — Vol. 197, I. 18. — P. 1-8.
44. Enzo-P: Petascale Enzo and the Cello Project j URL:
http:jjmngrid.ucsd.edujprojectsjcelloj (дата обращения: 25.07.2013).
45. PetaART: Toward Petascale Cosmological Simulations Using the Adaptive Refinement Tree (ART) Code j URL: http:jjwww.cs.iit.eduj zlanjpetaart.html (дата обращения: 25.07.2013).
46. Ferrari, A. A New Parallel SPH Method for 3D Free Surface Flows j A. Ferrari, et al. jj High performance computing on vector systems 2009. — 2010. — Part 4. — P. 179-188.
47. Van Straalen, B. Scalability Challenges for Massively Parallel AMR Applications j B. Van Straalen, J. Shalf, T. Ligocki, N. Keen, W. Yang jj In IPDPS ’09: Proceedings of the 2009 IEEE International Symposium on Parallel к Distributed Processing, Washington, DC, USA.
— P. 1--12.
48. Vshivkov, V. A Operator Approach for Numerical Simulation of Self-gravitation Gasdynamic Problem j V. Vshivkov, G. Lazareva, I. Kulikov jj Computational Technologies. — 200б. — Vol. 11, № 3. — P. 27-35.
49. Vshivkov, V. A Modified Fluids-in-cell Method for Problems of Gravitational Gas Dynamics / V. Vshivkov, G. Lazareva, I. Kulikov // Optoelectronics, Instrumentation and Data Processing. — 2007. — Vol. 43, I. 6. — P. 530-537.
50. Vshivkov, V. Computational Methods for Ill-posed Problems of Gravitational Gasodynamics / V. Vshivkov, G. Lazareva, A. Snytnikov, I. Kulikov, A. Tutukov // Journal of Inverse and Ill-posed Problems. — 2011. — Vol. 19, I. 1. — P. 151-166.
51. Aksenov, A.V. Symmetries and Relations Between Solutions of a Class of Euler-Poisson-Darboux Equations / A.V. Aksenov // Reports of RAS. — 2001. — Vol. 381, I. 2. — P. 176-179.
52. Vshivkov, V. Supercomputer Simulation of an Astrophysical Object Collapse by the Fluids-in-Cell Method / V. Vshivkov, G. Lazareva, A. Snytnikov, I. Kulikov // Lecture Notes in Computational Science. — 2009. — Vol. 5698. — P.414-422.
53. Petrov, M.I. Simulation of the Gravitational Collapse and Fragmentation of Rotating Molecular Clouds / M.I. Petrov, P.P. Berczik // Astronomische Nachrichten. — 2005. — Vol. 326. — P. 505-513.
Куликов Игорь Михайлович, к.ф.-м.н., младший научный сотрудник, Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук (Новосибирск, Российская Федерация), [email protected].
Черных Игорь Геннадьевич, к.ф.-м.н., научный сотрудник, Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук (Новосибирск, Российская Федерация), [email protected].
Глинский Борис Михайлович, д.т.н., заведующий лабораторией, Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук (Новосибирск, Российская Федерация), [email protected].
AstroPhi: A HYDRODYNAMICAL CODE FOR COMPLEX MODELLING OF ASTROPHYSICAL OBJECTS DYNAMICS BY MEANS OF HYBRID ARCHITECTURE SUPERCOMPUTERS ON Intel Xenon Phi BASE
I.M. Kulikov, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation),
I.G. Chernykh, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation),
B.M. Glinsky, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation)
There is a new hydrodinamical code AstroPhi for modeling of astrophysical objects dynamics on hybrid supercomputer proposed. This software package is optimized for using with Intel Xeon Phi calculations accelerators. AstroPhi code is based on combination of Godunov and author’s FlIC numerical methods for solving of gas dynamics equations. Fast Fourier Transform was used for Poisson equation solution. AstroPhi was tested on gas dynamics problems, Poisson equation solution and classical gravitational gas dynamics problems. The results of this tests and results
of gravitational collapse of astrophysical objects modeling proposed. The results of AstroPhi scalability based on Intel Xeon Phi runs are shown.
Keywords: numerical simulation, parallel computing, Intel Xeon Phi accelerated,
computational astrophysics.
References
1. Ardeljan N.V., Bisnovatyi-Kogan G.S., Moiseenko S.G. An Implicit Lagrangian Code for the Treatment of Nonstationary Problems in Rotating Astrophysical Bodies // Astronomy & Astrophysics. 1996. Vol. 115. P. 573-594.
2. Tutukov A., Lazareva G., Kulikov I. Gas Dynamics of a Central Collision of Two Galaxies: Merger, Disruption, Passage, and the Formation of a New Galaxy // Astronomy Reports.
2011. Vol. 55, I. 9. P. 770-783.
3. Gingold R.A., Monaghan J.J., Smoothed Particle Hydrodynamics: Theory and Application to Non-spherical Stars // Monthly Notices of the Royal Astronomical Society. 1977. Vol. 181. P. 375-389.
4. Luci L.B. A Numerical Approach to the Testing of the Fission Hypothesis // The Astrophysical Journal. 1977. Vol. 82, № 12. P. 1013-1024.
5. Collela P., Woodward P.R. The Piecewise Parabolic Method (PPM) for Gas-dynamical Simulations // Journal of Computational Physics. 1984. Vol. 54. P. 174-201.
6. Norman M. The Impact of AMR in Numerical Astrophysics and Cosmology // Lecture Notes in Computational Science and Engineering. 2005. Vol. 41. P. 413-430.
7. Hockney R.W., Eastwood J.W. Computer Simulation Using Particles New York: McGraw-Hill, 1981. 540 p.
8. Couchman H., Pearce F., Thomas P. Hydra Code Release / H. Couchman, F. Pearce, P. Thomas URL: http://arxiv.org/abs/astro-ph/9603116 (accessed: 25.07.2013).
9. Barnes J., Hut P. A Hierarchical O(Nlog N) Force-calculation Algorithm // Nature. 1986. Vol. 324. P. 446-449.
10. Dubinski J., Kim J., Park C., Humble R. GOTPM: A Parallel Hybrid Particle-Mesh Treecode // New Astronomy. 2004. Vol. 9. P. 111-126.
11. Fedorenko R. A Relaxation Method for Solving Elliptic Difference Equations // U.S.S.R. Computational Mathematics and Mathematical Physics. 1961. Vol. 1. P. 1092-1096.
12. Godunov S.K. A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations // Matematicheskii Sbornik. 1959. Vol. 47. P. 271-306.
13. Kulikovskii A.G., Pogorelov N.V., Semenov A.Yu. Mathematical Aspects of Numerical Solution of Hyperbolic Systems Moscow:Fizmatlit, 2001. 608 p.
14. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics Heidelberg:Springer-Verlag, 1999. 686 p.
15. Courant R., Isaacson E., Rees M. On the Solution of Nonlinear Hyperbolic Differential Equations by Finite Difference // Communications on Pure and Applied Mathematics. 1952. Vol. 5, № 3. P. 243-255.
16. Roe P. Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes // Journal of Computational Physics. 1997. Vol. 135, I. 2. P. 250-258.
17. Engquist B., Osher S.J. One-sided Difference Approximations for Nonlinear Conservation Laws jj Mathematics of Computation. 1981. Vol. Зб, № 154. P. 321-351.
18. Harten A., Lax P.D., Van Leer B. On Upstream Differencing and Godunov-type Schemes for Hyperbolic Conservation Laws jj Society for Industrial and Applied Mathematics Review. 1983. Vol. 25, № 1. P. 35-б1.
19. Einfeld B. On Godunov-type Methods for Gas Dynamics jj Society for Industrial and Applied Mathematics Journal on Numerical Analysis. 1988. Vol. 25, № 2. P. 294-318.
20. Batten P., Clarke N., Lambert C., Causon D.M. On the Choice of Savespeeds for the HLLC Riemann Solver jj Society for Industrial and Applied Mathematics Journal on Computing. 1997. Vol. 18, № б. P. 1553-1570.
21. Van Leer B. Towards the Ultimate Conservative Difference Scheme. V - A Second-order Sequel to Godunov’s Method jj Journal of Computational Physics. 1979. Vol. 32. P. 101-13б.
22. Wadsley J.W., Stadel J., Quinn T. Gasoline: a Flexible, Parallel Implementation of TreeSPH jj New Astronomy. 2004. Vol. 9, I. 2. P. 137-158.
23. Matthias S. GRAPESPH: Cosmological Smoothed Particle Hydrodynamics Simulations with the Special-purpose Hardware GRAPE jj Monthly Notices of the Royal Astronomical Society. 199б. Vol. 278, I. 4. P. 1005-1017.
24. Springel V. The Cosmological Simulation Code GADGET-2 jj Monthly Notices of the Royal Astronomical Society. 2005. Vol. Зб4, I. 4. P. 1105-1134.
25. Jin S., Xin Z. The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions jj Communications on Pure and Applied Mathematics. 1995. Vol. 48. P. 235-27б.
26. Godunov S.K., Manuzina Yu.D., Nazareva M.A. Experimental Analysis of Convergence of the Numerical Solution to a Generalized Solution in Fluid Dynamics jj Computational Mathematics and Mathematical Physics. 2011. Vol. 51. P. 88-95.
27. Ziegler U. Self-gravitational Adaptive Mesh Magnetohydrodynamics with the NIRVANA Code jj Astronomy к Astrophysics. 2005. Vol. 435. P. 385-395.
28. Mignone A., Plewa T., Bodo G. The Piecewise Parabolic Method for Multidimensional Relativistic Fluid Dynamics jj The Astrophysical Journal. 2005. Vol. 1б0. P. 199-219.
29. Hayes J. et al. Simulating Radiating and Magnetized Flows in Multiple Dimensions with ZEUS-MP jj The Astrophysical Journal Supplement Series. 200б. Vol. 1б5. P. 188-228.
30. Teyssier R. Cosmological Hydrodynamics with Adaptive Mesh Refinement. A New High Resolution Code Called RAMSES jj Astronomy к Astrophysics. 2002. Vol. 385. P. 337-3б4.
31. Kravtsov A., Klypin A., Hoffman Y. Constrained Simulations of the Real Universe. II. Observational Signatures of Intergalactic Gas in the Local Supercluster Region jj The Astrophysical Journal. 2002. Vol. 571. P. 5б3-575.
32. Stone J. et al. Athena: A New Code for Astrophysical MHD jj The Astrophysical Journal Supplement Series. 2008. Vol. 178. P. 137-177.
33. Brandenburg A., Dobler W. Hydromagnetic Turbulence in Computer Simulations jj Computer Physics Communications. 2002. Vol. 147. P. 471-475.
34. Schive H., Tsai Y., Chiueh T. GAMER: a GPU-accelerated Adaptive-Mesh-Refinement Code for Astrophysics jj The Astrophysical Journal. 2010. Vol. 18б. P. 457-484.
35. Murphy J., Burrows A. BETHE-Hydro: An Arbitrary Lagrangian-Eulerian Multidimensional Hydrodynamics Code for Astrophysical Simulations jj The Astrophysical Journal Supplement Series. 2008. Vol. 179. P. 209-241.
36. Springel V. E Pur Si Muove: Galilean-invariant Cosmological Hydrodynamical Simulations on a Moving Mesh jj Monthly Notices of the Royal Astronomical Society. 2010. Vol. 401. P. 791-851.
37. Bruenn S. et al. 2D and 3D Core-collapse Supernovae Simulation Results Obtained with the CHIMERA Code jj Journal of Physics. 2009. Vol. 180. P. 1-5.
38. Vshivkov V., Lazareva G., Snytnikov A., Kulikov I., Tutukov A. Hydrodynamical Code for Numerical Simulation of the Gas Components of Colliding Galaxies jj The Astrophysical Journal Supplement Series. 2011. Vol. 194, I. 47. P. 1-12.
39. Gonzalez M., Audit E., Huynh P. HERACLES: a Three-dimensional Radiation Hydrodynamics Code jj Astronomy к Astrophysics. 2007. Vol. 4б4. P. 429-435.
40. Krumholz M.R., Klein R.I., McKee C.F., Bolstad J. Radiation-hydrodynamic Simulations of the Formation of Orion-like Star Clusters. I. Implications for the Origin of the Initial Mass Function jj The Astrophysical Journal. 2007. Vol. бб7, I. 74. P. 1-1б.
41. Mignone A. et al. PLUTO: A Numerical Code for Computational Astrophysics jj The Astrophysical Journal Supplement Series. 2007. Vol. 170. P. 228-242.
42. Almgren A. et al. CASTRO: A New Compressible Astrophysical Solver. I. Hydrodynamics and Self-gravity jj The Astrophysical Journal. 2010. Vol. 715. P. 1221-1238.
43. Feng Y. et al. Terapixel Imaging of Cosmological Simulations jj The Astrophysical Journal Supplement Series. 2011. Vol. 197, I. 18. P. 1-8.
44. Enzo-P: Petascale Enzo and the Cello Project j URL:
http:jjmngrid.ucsd.edujprojectsjcelloj (accessed: 25.07.2013).
45. PetaART: Toward Petascale Cosmological Simulations Using the Adaptive Refinement Tree (ART) Code j URL: http:jjwww.cs.iit.eduj zlanjpetaart.html (accessed: 25.07.2013).
46. Ferrari A., et al. A New Parallel SPH Method for 3D Free Surface Flows jj High performance computing on vector systems 2009. 2010. Part 4. P. 179-188.
47. Van Straalen B., Shalf J., Ligocki T., Keen N., Yang W. Scalability Challenges for Massively Parallel AMR Applications jj In IPDPS ’09: Proceedings of the 2009 IEEE International Symposium on Parallel к Distributed Processing, Washington, DC, USA. P. 1-12.
48. Vshivkov V., Lazareva G., Kulikov I. A Operator Approach for Numerical Simulation of Self-gravitation Gasdynamic Problem jj Computational Technologies. 200б. Vol. 11, № 3. P. 27-35.
49. Vshivkov V., Lazareva G., Kulikov I. A Modified Fluids-in-cell Method for Problems of Gravitational Gas Dynamics jj Optoelectronics, Instrumentation and Data Processing. 2007. Vol. 43, I. б. P. 530-537.
50. Vshivkov V., Lazareva G., Snytnikov A., Kulikov I., Tutukov A. Computational Methods for Ill-posed Problems of Gravitational Gasodynamics jj Journal of Inverse and Ill-posed Problems. 2011. Vol. 19, I. 1. P. 151-1бб.
51. Aksenov A.V. Symmetries and Relations Between Solutions of a Class of Euler-Poisson-Darboux Equations jj Reports of RAS. 2001. Vol. 381, I. 2. P. 17б-179.
52. Vshivkov V., Lazareva G., Snytnikov A., Kulikov I. Supercomputer Simulation of an Astrophysical Object Collapse by the Fluids-in-Cell Method jj Lecture Notes in Computational Science. 2009. Vol. 5б98. P.414-422.
53. Petrov M.I., Berczik P.P. Simulation of the Gravitational Collapse and Fragmentation of Rotating Molecular Clouds jj Astronomische Nachrichten. 2005. Vol. 32б. P. 505-513.
Поступила в редакцию їв октября 2Q13 г.