CC BY
15ASTROPHI 2.0: НОВЫЙ КОД ВЫСОКОГО ПОРЯДКА
ТОЧНОСТИ ДЛЯ ГИДРОДИНАМИЧЕСКОГО МОДЕЛИРОВАНИЯ АСТРОФИЗИЧЕСКИХ ТЕЧЕНИЙ НА ГИБРИДНЫХ СУПЕРЭВМ, ОСНАЩЕННЫХ УСКОРИТЕЛЯМИ INTEL XEON PHI
И. M. Куликов, И. Г. Черных
Институт вычислительной математики и математической геофизики СО РАН,
630090, Новосибирск, Россия
УДК 519.6, 524.3
В статье представлен новый гидродинамический программный код AstroPhi 2.0 для численного моделирования астрофизических процессов на гибридных суперЭВМ, оснащенных ускорителями Intel Xeon Phi. Описаны детали параллельной реализации кода и элементы со-дизайна численного алгоритма, которые позволили сделать эффективную программную реализацию. В рамках одного ускорителя было получено 134-кратное ускорение, 92-процентная эффективность была получена при использовании 64 ускорителей. С помощью данного кода была смоделирована задача столкновения галактик.
Ключевые слова: математическое моделирование, суперкомпьютерные вычисления, параллельные вычислительные методы, вычислительная астрофизика, столкновение галактик, ускорители Intel Xeon Phi.
A new hydrodvnamic code AstroPhi for numerical simulation an astrophvsical processes on hybrid supercomputers by means Intel Xeon Phi accelerator 2.0 was presented. The details of parallel development and efficiency co-design elements of numerical algorithm was described. A speed-up of 134 times was obtained within a single Intel Xeon Phi accelerator. The use of 64 Intel Xeon Phi accelerators resulted in 92 % parallel efficiency. By means AstroPhi 2.0 code a problem of interacting galaxies was simulated.
Key words: numerical modeling, parallel computing, parallel numerical method, computational astrophysics, interacting galaxies, Intel Xeon Phi accelerators.
Введение. Предметом изучения современной астрофизики является исследование физических процессов во Вселенной, как эти процессы повлияли на самоорганизацию и эволюцию астрофизических объектов, а также на дальнейшую их динамику и взаимодействие. Существенность учета гравитационного и магнитного полей, а также сложность воспроизведения космических условий в лабораторных условиях накладывают значительные ограничения на экспериментальное изучение астрофизических объектов. По существу, математическое моделирование — это основной, если не единственный, подход к теоретическому исследованию астрофизических процессов и объектов.
В контексте космологического моделирования скоплений галактик в рамках рассматриваемой стандартной космологической модели ЛСБМ необходимо моделирование прежде
Работа поддержана грантами Российского фонда фундаментальных исследований № +15-31-20150 мол-а-вед, № 15-01-00508 и № 14-01-31199, грантом Президента РФ № МК-6648.2015.9.
всего холодной темной материи, для описания которой используется бесстолкновительная модель N-тел, За последние годы в рамках этой модели был проведен ряд глобальных космологических моделирований [1, 2, 3]. Однако, для перехода к следующему уровню, к отдельным галактикам, такой подход имеет значительные ограничения из-за невозможности получения ряда эллиптических галактик. Поэтому необходимо моделировать также барионную материю в гидродинамической модели и проводить совместное космологическое моделирование, такое как проекты „Bolshoi" [4], „Illustris" [5, 6, 7] и „Eagle" [8, 9], Уточнение космологической модели происходит не только в сторону использования многофазной модели, но и в сторону уеложения так называемой subgrid (подсеточной) физической модели.
Структурным элементом скопления галактик являются собственно сами галактики. Задачи моделирования динамики галактик можно разделить по времени их динамики. Так эволюция отдельной галактики составляет до нескольких миллиардов лет, в то время как взаимодействие отдельных галактик составляет несколько сотен миллионов лет. Основными проблемами являются вопрос о профиле плотности в галактиках, скорости процесса звездообразования, а также скорости вращения галактик. Проблема профиля скорости вращения поднимается в связи с анализом отношения Tully-Fisher и распределением темной материи в галактиках и гравитационного потенциала [10]. В исследовании столкновения галактик отдельно стоит выделить проект GALMER парижской обсерватории — базу данных вычислительных экспериментов [11] по столкновению различных типов галактик, В задачах столкновения галактик [12] процессы (звездообразование [13], AGN [14], образование сверхмассивных двойных черных дыр [15, 16], химокинетика [17]) значительно ускоряются, и необходим явный учет этих процессов в математической модели в виде уже озвученной ранее subgrid (подсеточной) физической модели.
Все больший интерес представляют динамики молекулярных облаков, их химическая эволюция [18, 19] и химическая эволюция карликовых галактик [20, 21, 22], которые могут быть сформулированы в одном контексте в связи со сравнимыми линейными размерами и массами. Так была смоделирована химическая эволюция дисковых галактик в зависимости от различного распределения массы [23], в эволюции Млечного Пути [24, 25], в исследовании металличноети дисковых галактик [26], в эволюции галактик с галактическим ветром [27], в галактическом гало [28], Особый интерес представляют химическая эволюция молекулярных облаков и образование сложных элементов вплоть до спиртов [29].
Основным фокусом наших исследований является задача столкновения галактик. Относительно особенностей этой задачи и строилась триада „математическая модель — чис-
"
мический код AstroPhi 2.0. Одной из главных проблем моделирования взаимодействующих галактик является соотношение масштабов. Так масса одной галактики составляет 1013M©, размер 10 килопарсек, а масса рядовой звезды — порядка M© и размер 10-7 парсек, что приводит к разрыву в 13 порядков для массы и 14 порядков для размера. Данное обстоятельство приводит к необходимости использования мощнейших из доступных суперЭВМ. На сегодняшний день самые производительные суперкомпьютеры построены на гибридной архитектуре на основе графических ускорителей и ускорителей Intel Xeon Phi. Так в ноябрьской версии 2014 года Тор500 первые два суперкомпьютера (и половина ТорЮ) построены на этих технологиях. Очевидно, что на основе именно гибридных архитектур будет построен первый экзафлопеный суперкомпьютер. Конечно, гибридные
Таблица 1
Достоинства и недостатки графических ускорителей и ускорителей Intel Xeon Phi
Достоинства NVIDIA Tesla/Kepler
Достоинства Intel Xeon Phi
Большое число ядер Методическая поддержка Контроль над оптимизацией Контроль передачи данных Большое число библиотек
Простота использования Стандарт ОрепМР Реализация операции редукции
Переносимость Поддержка стандарта IEEE 754
Недостатки NVIDIA Tesla/Kepler
Недостатки Intel Xeon Phi
Сложность программирования Ограничения на алгоритмы Ограничения на структуры данных Отсутствие редукции Непереносимость кода Потеря точности вычислений Слабая производительность ядер
Малое число ядер Слабая методическая поддержка Слабая производительность ядер
архитектуры имеют как достоинства, так и недостатки, которые мы попытались перечислить в табл. 1,
Общим устранимым недостатком является слабая производительность ядер, что заставляет достаточно эффективно использовать именно независимый способ организации вычислений. Особенности архитектуры, сетевой инфраструктуры, организации памяти и вычислений делают разработку программного обеспечения для таких суперЭВМ сложной научной задачей, которая требует детальной проработки от физической постановки задачи до инструментов разработки, В этом случае мы говорим о концепции со-дизайна высокопроизводительных вычислений [30], В связи с этим модель взаимодействующих галактик была сформулирована в виде системы гиперболических уравнений как для описания газовой компоненты, так и для описания бесстолкновительной компоненты, описывающей звезды и темную материю [31], Такая формулировка позволяет нам использовать согласованный подход к фазовому переходу необходимого для моделирования процесса звездообразования и взрыва сверхновых звезд, единый численный метод высокого порядка точности, способ задания структур данных и параллельного алгоритма,
В последние два десятилетия из широкого диапазона газодинамических численных методов для решения нестационарных трехмерных астрофизических задач используются два основных подхода. Это лагранжев подход, в основном представленный SPH-методом [32, 33] (Smoothed Particle Hydrodynamics), и эйлеров подход с использованием адаптивных сеток или АМН [34] (Adaptive Mesh Refinement). В последние пять лет появился ряд программных пакетов с использованием комбинации лагранжева и эйлерова подходов. Основной проблемой SPH-метода являются поиск соседей и организация их гравитационного взаимодействия. Для эффективного решения этой задачи был разработан ряд алгоритмов, Например, particle-particle/particle-mesh или Р3М метод [35], адаптация Р3М метода с использованием иерархичности расчетной сетки АР3М [36], tree алгоритм [37], комбинация tree алгоритма и particle-mesh подхода (Тгее-РМ метод) [38], Для решения уравнения Пуассона в сеточных методах используются в основном метод сопряженных градиентов (CGM), метод быстрого преобразования Фурье (FFT), метод последовательной верхней
релаксации (SOR) и метод Федореико [39] или многосеточный метод (MG).
Для численного решения газодинамических задач широкое применение получил метод Годунова [40], основным структурным элементом которого является задача о распаде произвольного разрыва (задача Римана) с параметрами газа в соседних ячейках разностной сетки. Как правило, параметры газа в соседних ячейках достаточно близки, что создает благоприятные условия для применения упрощенного алгоритма решения задачи о распаде разрыва. Различные алгоритмы получения приближенного решения задачи Римана дали большой класс методов [41, 42]. Основными методами являются методы типа Куранта-Изаксона-Риса [43] и Роу [44], которые строятся на основе использования различным образом линеаризованных гиперболических систем уравнений, Ошера [45], где решение задачи Римана строится только с использованием волн Римана, Основоположным подходом к оценке скоростей волн является двух вол новой метод Хартепа-Лакса-Вап Лира (известный в литературе как HLL) [46], в котором учитываются левые и правые разрывы без рассмотрения контактного разрыва, В схеме HLL, использующей консервативные переменные, предложен простой, но эффективный способ выбора скоростей движения этих волн по максимальным наклонам характеристик в соседних ячейках разностной сетки. При этом веер волн разрежения заменяется скачком, но со скоростью распространения,
соответствующей максимальному наклону характеристик в этой волне разрежения, По"
теристик. Расчет ударных волн также проводится со скоростью, превышающей точное значение, В этой связи схема HLL эффективна при расчете ударных волн и зон разрежения, Однако в расчете энтропийных скачков методом установления принятое допущение
"
дификации HLL, такие как HLLE [47], где особо учитываются крайние собственные числа линеаризованной задачи, и HLLC [48], где производится дополнительный учет центрального разрыва, движущегося со скоростью, равной центральному собственному значению линеаризованной задачи Римана,
Также на основе метода Годунова были сделаны реализации высокого порядка — монотонная противопотоковая схема второго порядка точности Ml SCI. [49] и TVD схемы [50], третьего порядка кусочно-параболический метод РРМ [51], неосциллпрующий схемы с весами WENO [52], Однако, что понимается под высоким порядком точности в случае разрывных решений, не очень понятно [53], Основной недостаток методов высокого порядка — это необходимость введения ограничителя потока (flux limiter в зарубежной литературе). Такой подход зачастую приводит к искажениям численного решения, например, к неустранимой ошибке скорости ударной волны. Также существенным недостатком методов высокого порядка точности является необходимость использования отличного от компактного сеточного шаблона вычислений, что приводит к необходимости использования большего числа слоев перекрытия, что приводит к увеличению нагрузки на сетевую инфраструктуру и, как следствие, падению производительности, В связи с этим, для увеличения порядка точности нашего метода мы использовали модификацию метода РРМ — метод Pieeewise-Parabolie-Method on local stencil (PPML) [54, 55],
В рамках лагранжева подхода на основе SPH метода были разработаны пакеты Hydra [56], Gasoline [57], GrapeSPH [58], GADGET [59], В рамках эйлерова подхода (в том числе и с использованием AMR) были разработаны пакеты NIRVANA [60], FLASH [61], ZEUS [62], ENZO [34], RAMSES [63], ART [64], Athena [65], Pencil Code [66], Heracles [67], Orion [68], Pluto [69], CASTRO [70]. Эйлеров подход с использованием AMR был впервые ис-
Таблица, 2
Достоинства и недостатки SPH и AMI! методов
Достоинства БРН: Достоинства AMI!:
Робастность алгоритма Инвариантность решения Простота реализации Адаптация под геометрию задачи Высокая точность потенциала Апробированная методология Отсутствие искусственной вязкости Корректность на ударной волне Воспроизведение разрывов Воспроизведение неустойчивостей Воспроизведение турбулентности
Недостатки БРН: Недостатки AMI!:
Искусственная вязкость Радиус сглаживания Осцилляции на разрывах Подавление неустойчивостей Слабая масштабируемость Сложность реализации Сеточные эффекты Проблема разрешения сетки Проблема инвариантности Слабая масштабируемость
пользован на гибридных суперкомпьютерах, оснащенных графическими ускорителями, в пакете GAMER [71]. В табл. 2 приведены достоинства и недостатки этих методов.
Пакеты BETHE-Hvdro [72], AREPO [73], CHIMERA [74], GIZMO [75] и авторская реализация PEGAS/GPUPEGAS/AstroPhi [76, 77, 78] основаны на комбинации лагранжева и эйлерова подходов и в различной мере устраняют недостатки при сохранении достоинств базовых методов. В табл. 3 приведены методы и технологии, использованные для реализации основных программных пакетов.
В работе описаны математическая модель взаимодействующих галактик с учетом различных subgrid-physics процессов и численный метод, используемый для ее решения. Приведены схема параллельной реализации и анализ производительности различных подсистем кода AstroPhi 2.0, а также продемонстрированы результаты вычислительных экспериментов по столкновению галактик. Несмотря на то, что основная область применения кода — это задачи столкновения галактик, код может с успехом использоваться и в других астрофизических приложениях, где необходимо гидродинамическое моделирование.
1. Описание математической модели. Основными ингредиентами галактик являются столкновительная компонента, моделирующая газовую составляющую галактик и равномерно распределенную в ней пыль, и бееетолкновительная компонента, моделирующая звездную компоненту и темную материю. Для учета химических реакций будем рассматривать уравнения многокомпонентной односкоростной газовой динамики. Для описания бееетолкновительной компоненты будем использовать уравнения для первых моментов бееетолкновительного уравнения Больцмана. Такой подход был успешно апробирован на моделировании эволюции [31] и столкновении галактик [77], а также он позволяет термодинамически согласованно описать процессы фазового перехода между компонентнами, которые происходят при звездообразовании и взрыве сверхновых.
1.1. Уравнения многокомпонентной многофазной гидродинамики. Для описания газовой компоненты будем использовать систему уравнений односкоростной компонентной гравитационной газовой динамики, записанную в эйлеровых координатах:
-£ + V • (pu) = S-D,
Таблица 3
Основные характеристики астрофизических пакетов
Код Численный метод Технологии
Hydra SPH + Adaptive P3M + FFT HPF
Gasoline SPH + Tree code + Multipole Method MPI
GrapeSPH SPH + Direct Summation GRAPE
GADGET-2 SPH + TreePM + FFT MPI
NIRVANA AMR+HLL + Multigrid MPI
FLASH AMR+PPM + Multigrid MPI
ZEUS-MP Finite difference method + FFT+Multigrid MPI
ENZO AMR+PPM + FFT+Multigrid MPI
RAMSES AMR+HLLC + Multigrid+CG OpenMP+MPI
ART AMR+MUSCL + FFT MPI
Athena Roe's solver + FFT MPI
Pencil Code Finite difference method + FFT HPF+MPI
Heracles MUSCL + CG MPI
Orion AMR+MUSCL + Multigrid _
Pluto AMR+HLLC + Analytical MPI
CASTRO AMR+PPM + Multigrid MPI+OpenMP
GAMER AMR+TVD + FFT+SOR CUDA+MPI
BETHE-Hvdro Arbitrary Lagrangian-Eulerian + Matrix Inverse _
AREPO Moving mesh + MUSCL + TreePM + FFT MPI
CHIMERA Moving mesh + PPM + Analytical _
PEG AS FlIC+Godunov+PPML + FFT MPI+CUDA+OpenMP
дРн , ^ сря „Рн
-тгг + V ■ (рни) = -Sh,h2 + S--D—,
dt p p
ддН + V- (Рн U) = SH,H2 + SРя2 -DРя2 dt p p
dpu ~dt~
+ V ■ (puu) = -Vp - рУ(Ф) + vS - uD,
^ + V ■ (pEu) = -V ■ (pv) - (рУ(Ф),и) - Л + Г - eD, dt p
dpe D
+ V ■ (peu) = -(7 - 1WV ■ и - Л + Г - e—, dt p
12
pe = 2 Pu + Pe,
P = (Y - 1)pe.
Для описания бесстолкповительпой компоненты будем использовать систему уравнений для первых моментов бесстолкновительного уравнения Больцмана, записанную также в эйлеровых координатах:
dt + V-(nv) = D-S,
дплт
-— + V ■ (пуу) = -УП - пУ(Ф) + и— - Ув,
дг
дрШ V
— + V ■ (рШ?) = -V ■ (Пу) - (пУ(Ф),У) - Г + е-, сл р
ддт + У(П« у) = -2ПУи - Г+е -
Ттг 1 ->2 : Пхх + Пуу + П.г.г
рш = 2 ру2 +-^-,
уравнение Пуассона для обеих компонент записывается в виде:
ДФ = 4пС(р + п),
где р — давление газа, рН — плотность атомарного водорода, рН2 _ плотность молекулярного водорода, 8Н,н2 — скорость образования молекулярного водорода из атомарного, р = рН+рН2 _ плотность смеси газа, п — плотность бесстолкновительной компоненты, и — скорость газовой компоненты, V — скорость бесстолкновительной компоненты, рЕ — плотность полной механической энергии газа, рШ — плотность полной механической энергии бесстолкновительной компоненты, Ф — гравитационный потенциал, е — плотность внутренней энергии газа, 7 — эффективный показатель адиабаты, П^ = (Пхх, Пуу, Пхх) — диагональный тензор дисперсии скоростей бесстолкновительной компоненты, в — скорость образования сверхновых звезд, — — скорость звездообразования, Л — функция компто-
Г
1,2, Модель „подсеточной" физики. Мы остановимся на следующих, на наш взгляд, основных „подееточных" процессах, играющих основную роль при столкновении галактик:
1) самосогласованная модель химической кинетики, В нашей модели мы остановимся на процессе образования молекулярного водорода;
2) процесс звездообразования, взрыва молодых звезд и связанного с последним процесса нагревания;
3) комптоновское охлаждение газа.
Остановимся на каждой модели подробнее,
1,2,1, Образование молекулярного водорода. Молекулы водорода в межгалактическом пространстве формируются на поверхности частиц и диссоциируют космическим излучением, Предполагая, что плотность газа пропорциональна плотности частиц из-за хорошего перемешивания частиц и газа в нашей модели, концентрация молекулярного водорода определяется следующим выражением [79]:
¿пН2
¿г
Ядт (Т )пн (Пн + 2ПН2 ) - (Сн + &гзз^Н2 ,Ау )) ПН2
где пн2 и пН — концентрация молекулярного и атомарного водорода, МН2 — столбцевая плотность молекулярного водорода, скорость образования молекулярного водорода на пыли задается функцией [80]:
Ядг(Т) = 2,2 х 10-18Б^-1
где Б = 0,3 — эффективность образования молекулярного водорода на пыли [81], скорость ионизации водорода за счет космических лучей задается функцией [82, 83]:
<н = 6 х 10-18в-1,
где Ль — погашение [84], Скорость фотодиееоциации <dгss(Nн2,Ау) записывается в виде [85]:
,лу ) =
(0)/вhieИ(N (Н^/^Лу ),
где
<&5в(0) = 3,3 х 1,7 х 10-115-1 — неэкранированная скорость фотодиссоциации [86],
) = ехр(-т^1000(Лу)) — скорость абеорпции па пыли [85], где т^,1000(Лу) = 3,74Лу = 10-21 (М(Н) + N(Н2)) — оптическая глубина на частицах пыли на длине волны Л = 1000 А, где N(Н) и N(Н2) — столбцевые плотности. Функцию коэффициента самозащиты можно аппроксимировать [85]:
/выем^(Н2)) = п ",9,, ,2 + 7°=ехр (-8,5 х 10-4^Г+Х) , (1+ х/Ь.5)2 л/1 + х V /
0,965 0,035 (1 + х/65)2 + где
х = N(Н2)/5 х 1010т2,
Ь5 = Ь/107 м/с, где Ь — параметр Допплеровекого расширения. После вычисления концентрации молекулярного водорода пн2, определяется эффективный показатель адиабаты
2
5п# + 7п#2
7
3п# + 5п#2
Уравнение для описания изменения концентрации водорода имеет аналитическое решение, которое и было использовано в реализации кода Ав^оРЫ 2,0,
1,2,2, Процесс звездообразования и образование сверхновых звезд. Будем использовать следующее необходимое условие процесса звездообразования, сформулированное в виде [87]:"
Т < 104К, у ■ и < 0, р > 1,64
М®
рс-3
тогда скорость звездообразования может быть сформулирована в виде:
^ «п „ р „ 3/2 /32С
® = - ==С р3/У ат
тауп V 3п
где С = 0,034 — коэффициент эффективности звездообразования. Скорость образования сверхновых звезд записывается в виде [88]:
5 = Д™ = «Р = вС— = вС п3/2./
т¿у„ V 3п
где в = 0,1 — коэффициент взрыва молодых звезд. При взрыве одной звезды солнечной
1051
сверхновых может быть записана в виде:
51 MSN Г= 1051 — erg, MQ g
MSN
1,2,3, Комптоновекая функциях охлаждения. Галактический газ, разогретый за счет столкновения до температуры - 104 - 108K, описывается функцией охлаждения [89]:
Л ~ 10-22nHcm-3erg.
где nH — концентрация атомарного водорода,
2. Численный метод решения. Для решения гидродинамических уравнений был использован оригинальный численный метод, основанный на комбинации operator splitting подхода, метода Годунова, схемы Рое и кусочно-параболической реконструкции на локальном шаблоне. Такой метод объединяет все достоинства перечисленных методов и обладает высокой степенью параллелизации, на чем остановимся подробнее в следующем разделе. Для решения уравнения Пуассона используется метод, основанный на быстром преобразовании Фурье,
2,1, Эйлеров этап решения гидродинамических уравнений. На эйлеровом этапе численного метода решаются уравнения без учета адвективных членов и функций, моделирующих ,,подееточные" процессы. Таким образом, мы будем рассматривать уравнения для газовой компоненты в виде:
дри дрЕ дре _
~дГ = -vp ~дТ = (pv), иг = -(Y - 1)рбУ"и,
для бееетолкновительной компоненты в виде:
dnv = -vn> ddW = -V- (nv). ^ = --V- и.
Для аппроксимации пространственных производных используется решение линеаризованной задачи Римана, Для этого на всех интерфейсах между ячейками (L — обозначение левой ячейки, R — обозначение правой ячейки) происходит осреднение скорости и давления по формулам:
= рГ + рд 2 = Pl/Pl + Pr/Pr
/Рь + /Pr. /рь + VPr '
Такой способ осреднения связан с необходимостью более аккуратного учета границы газ-вакуум, чем то, как это сделано в оригинальной схеме Рое [44], При этом будем предполагать, что в рассматриваемых ячейках решение является кусочно-параболической функцией. Процедура построения локальных парабол подробно описана в приложении настоящей статьи, а ее анализ и детали построения описаны в классической работе [51], Эта процедура нам понадобится также в следующем этапе метода.
Опустив подробные выкладки для вывода схемы с помощью метода Годунова (численный метод и его верификация подробно описаны в работе [90]), решение задачи Римана для эйлерова этапа решения уравнений для газовой компоненты суть:
U = ul(-M)+ ur(At) + Pl(-Aí) - Pr(Xí) I (^J~pL + v/pR)S
2 2 V yl+YR (p3/2 + з/2
(PL + Pr )(Pl/Pl + Pr/Pr)
2
= PL(-At)+ pfl(At) + Ub(-At) - ufl(At) /^(pf + pR/2)(plVpZ + Рд VpR)
2 2 V (VpL + VpR)2
для бееетолкновительной компоненты суть:
V = ^¿(-^t) + + nL(-jut) - Пд(^) / (^nL + ^nR)2
2 2 V 3(ni/2 + 4/2)(nL^tL + Пд^ПД)
П = Щ(-^) + Пд(^t) + ^¿(-^t) - Vfl(^t) 3(n3j/2 + пД/2+ Пд^йД)
= 2 2 V (VnL + VnR)2 ,
где
Л = /Д+1Д (рьТР! + Ряур?) = /3(П£^Т + Пнуп?)
V рГ + рГ ' * V "х2 + 4/2 ' = 9? - | (Д* - «6 (1 - ^) ) • = ^ + 21 (Д* + «К1 - Ж
Приведем подробную процедуру построения параболы и параметров Д^, «6, Для
определенности будем конструировать кусочно-параболическую функцию произвольного параметра «(х) на регулярной сетке с шагом на интервале [х^-1/2,хг+1/2]. В общем виде парабола может быть записана как:
(6),
q(x) = qL + С (A® + q(6)(1 - С)) ,
где qi — значение в центре ячейки, С = (x-х-г^, Адг = qL - qR и q(6) = 6(q - 1/2(qL + qR)) ПрИ уСловии сохранения консервативности, то есть:
il
г xi+1/2
qi = h / q(x)dx.
xi-1/2
Для конструирования значений qR = qL+1 = qi+i/2 будем использовать интерполяционную функцию четвертого порядка точности:
qi+i/2 = 1/2(qi + qi+1) - 1/6(Sqi+1 - Sqi),
где Sqi = 1/2(qi+1 - qi-1). Далее опишем алгоритм получения локальной параболы. На вход алгоритма подаются значения в точках ячеек q^ На выходе алгоритма определяются все параметры кусочно-параболических функций на всех интервалах [xi-1/2,xi+1/2].
Шаг 1, На первом шаге мы конструируем значения Sqi = 1/2(qi+1 - qi-1). Для этого нам необходимо знание только соседних ячеек qi+1,qi - 1, Для избежания экстремумов
Sqi
{mm(|Sqj|, 2|qj+1 - qi|, 2|qj - qi-1 |)sign(Sqi), (qi+1 - qi)(qi - qi-1) > 0 0, (qi+1 - qi)(qi - qi-1) ^ 0
В случае параллельной реализации на архитектурах с распределенной памятью мы должны сделать обмены одного слоя перекрытия расчетной области средствами MPI, После чего пересчитываем значения на границе с помощью интерполянта четвертого порядка точности:
qR = qi+i = qi+1/2 = 1/2(qi + qi+i) - 1/6(5mqi+i - 5mqi).
Шаг 2, На втором шаге алгоритма мы начинаем конструировать саму локальную параболу с помощью формулы:
Aqi = qL - qR, q(6) =6(qi - 1/2(qL + qf)).
В случае немонотонности локальной параболы (такое имеет место на разрывах) мы перестраиваем значения на границах qL,qR по формулам:
qL = qi, qf = qi, (qf - qi)(qi - qf ^ 0,
qL = 3qi - 2qf, Аад(6) > (Aqi)2, qf = 3qi - 2qf, Аад(6) < -(Aqi)2.
Таким образом, граничные значения удовлетворяют условиям монотонности.
Шаг 3, На третьем шаге перестроим параметры параболы с учетом новых значений на границах ячеек:
Aqi = qL - qR, q(6) = 6(qi - 1/2(qf + qf)).
В результате локальная парабола в каждой ячейке [xi-1/2,xi+1/2] получена. Стоит отметить, что параболы могут иметь разрыв на границах ячеек, что в случае использования классического кусочно-параболического метода (РРМ) приводит к необходимости решения задачи Римана для парабол, В нашем случае локальные параболы используются как составная часть задачи Римана,
2.2, Лагранжев этап решения гидродинамических уравнений. На лагранжевом этапе происходит адвективный перенос гидродинамических параметров, и все уравнения на лагранжевом этапе имеют вид:
ж + v (fv) = 0-
где f может быть плотностью р, и, импульсом pu, nv, плотностью общей механической pE, nW или внутренней ре, П^ энергии газа. Для решения уравнений используется аналогичный, что и на эйлеровом этапе, подход. Для вычисления потока F = fv при А = |v| используется формула:
F = v Л ,
\ fR(At),v< 0
где fb(-At) и ff(At) — кусочно-параболические функции для величины f и скорость на интерфейсе между ячейками вычисляется по усреднению Рое:
= Vb^pL + Vf^pR
лfpl + VpR .
Для построения кусочно-параболического решения используется процедура, которая приведена в предыдущем разделе,
2.3, Учет „подсеточных" процессов. На третьем этапе численного метода решаются уравнения:
др дГ = = 5-—,
дрн дГ = - 5н,н2 + 5 РН -Р х>н Рн Р
дРя2 дГ = = 5н,н2 + 5 Рн Р — РН2 , Р
^ = -рУ(Ф) + - и—,
^ = -(рУ(Ф),и) - Л + Г - е-, дг р
дре д тл —
-1- = -Л + Г - е-, дГ р
дп ^
ж =—-5 ■
дпу
= -гсУ(Ф) + и— - ,
ддЕ = _(„у(Ф),*) - Г + е-,
дг р
дПг! = -Г + е -,
дГ 3р'
с помощью метода Эйлера для решения обыкновенных ОДУ,
2.4, Использование переопределенной системы, уравнений. На завершающем этапе решения гидродинамических уравнений происходит корректировка решения, В случае границы газ-вакуум — с использованием формулы [91]:
|у| = - б), (Е - V2/2)/Е ^ 10-3,
в остальной области происходит корректировка, которая гарантирует неубывание энтропии [92]:
|рб| = (р£ - ^ , (Е - ^2/2)/Е < 10-3.
Такая модификация обеспечивает детальный баланс энергий и гарантирует неубывание энтропии,
2.5, Решение уравнения Пуассона. После решения гидродинамических уравнений необходимо восстановить гравитационный потенциал по плотностям газовой и бееетолкнови-тельной компонент. Для этого мы будем использовать 27-точечный шаблон для аппроксимации уравнения Пуассона:
38 4
- уФгД! + 9 (Ф»-1,к,г + Фг+1Д1 + Ф»,к-1,г + Фг,к+1>г + ФгД1-1 + ФгДш) +
1 (Фг-1,й-1,г + Фг+1,й-1,г + Фг-1,й+1,г + Фг+1,й+1,г + Фг-1,й,г-1 + Фг+1,й,г-1 + У
Фг-1,к,1+1 + Фг+1,к,1+1 + Фг,к-1,1-1 + Фг,к+1,1-1 + Фг,к-1,1+1 + Ф
— (Фг-1,к-1,1-1 + Фг-1,к+1,1-1 + Фг-1,к-1,1+1 + Фг-1,к+1,1+1 + Фг+1,к-1,1-1 + Фг+1,к+1,1-1 +
36
Фг+1,й-1,г+1 + Фг+1,й+1,г+1) = 4пСН2 (рг,к,1 + щ,к,1).
Это связано с тем, чтобы обеспечить максимальную инвариантность решения относительно поворота. Алгоритм решения уравнения Пуассона будет состоять из нескольких этапов.
Этап 1, Постановка граничных условий для уравнения Пуассона, Для постановки граничных условий для гравитационного потенциала па границе области В будем использовать первые члены мультипольного разложения — статический, осевой и центробежный моменты инерции:
М М
Ф(х,у,г)Ь =----3 (4 + 1у + 12 - 31о),
г г з
где
I = (х21Х + у% + г2 4) - 2 (ху1ху + хг!хг + уг1уг)
Ъ = ^ (г2 + у2) т, 1у = ^ (х2 +3 т, Ъ = ^ (х2 + У2) т,
3 3 3
1ху 'У ] х3 у3 т3, ^ ] х3 г3 т3, 1уг ^ ] у3 23 т3,
333
где х, у, г — координаты центра ячеек на границе области, г = \]х2 + у2 + г2 — расстояние до центра области, х3, уу, 23 — координаты очередной ячейки, т3 — масса очередной М
лены, то они подставляются в 27-точечный шаблон для определения итоговой плотности рг,к,1 + пг,к,1 на границе В.
Этап 2, Преобразование плотности в пространство гармоник. Итоговая плотность представляется в виде суперпозиции по собственным функциям оператора Лапласа:
I 1пкт т!п
Рг,к,1 + Пг,к,1 = &3тпехр I —р + ^ +
3тп \
где I, К, Ь — число ячеек то каждой координате, 1 — мнимая единица. Для этого используется быстрое преобразование Фурье,
Этап 3, Решение уравнения Пуассона в пространстве гармоник. Мы предполагаем, что потенциал также представлен в виде суперпозиции по собственным функциям оператора Лапласа:
/ 1пкт т!п
Фг,к,1 = Ф3тп&хр\ —р + к +
3тп \
При подстановке такого разложения в схему аппроксимации уравнения Пуассона мы получаем достаточно простую формулу для вычисления амплитуд гармоник потенциала:
Лагранжев этап (80 %)
Рис. 1. Процентное; соотношение между этапами в коде AstroPhi 2.0
jmn 1 ( 2W( 2WЛ 2W^
После чего необходимо проделать обратное быстрое преобразование Фурье гармоник потенциала в функциональное пространство гармоник.
3. Описание параллельной реализации. Со-дизайн 1301 физико-математической модели, численного метода и структур данных позволяет использовать геометрическую декомпозицию расчетной области с одним слоем перекрытия подобластей. Такую возможность мы имеем за счет построения парабол па предыдущем шаге, что требует только локального взаимодействия между ячейками. На рис. 1 приведены процентные соотношения между этапами.
Дня решения уравнения Пуассона, в основе которого быстрое преобразование Фурье дня суперЭВМ с распределенной памятью, была использована библиотека FFTW |93|, В основе этой библиотеки лежит процедура ALLTOALL, которая „транспонирует" трехмерный массив, перераспределяя значительные объемы памяти между всеми процессами. Безусловно, это дорогая сетевая операция, которая требует отказа от всего алгоритма в случае использования сколь либо значительного количества вычислителей. Однако эта процедура в случае использования сетевой инфраструктуры InfiniBand пе занимает критическое время и, по всей видимости, оптимизирована па низком сетевом уровне |94|,
Основными этапами вычислительной схемы являются эйлеров и лагранжев этапы. Мы сосредоточимся имешю па этих этапах как па наиболее затратных. Из оставшихся этапов наибольшее время занимает процесс сохранения результатов па диск, что, в принципе, ожидаемо, и далее мы приведем исследования скорости записи па диск контрольных точек. Также вне нашего рассмотрения в плане ускорения останутся процедуры, в которых "
фактически происходит копирование памяти из одной области в другую, в дальнейшем мы также рассмотрим эти операции отдельно с точки зрения обобщенной функции MEMCPY.
Отдельно остановимся па процедуре вычисления шага по времени, исходя из условия Kvpairra, В случае использования графических ускорителей данная процедура была реализована только па CPU (также было сделано и в коде GAMER |71|). Причина этого -отсутствие эффективной реализации редуцирующей операцией min в технологии CUDA, в то время как в ОрепМР такая операция эффективно реализована. Стоимость этой процедуры составляет порядка одного процента от общего времени вычислений и практически
/ 7í"
MPI
#pragma omp parallel for
Рис. 2. Декомпозиция расчетной области в коде AstroPhi 2.0. Серым цветом обозначены слои
перекрытия подобластей
по влияет па эффективность параллельной реализации. Однако при увеличении количества графических ядер до нескольких тысяч и стократном ускорении в рамках одного графического процессора суммарно всех остальных процедур может возникнуть курьез-пая ситуация, когда процедура вычисления шага по времени будет выполняться дольше всех остальных. При том, что авторами уже было достигнуто 55-кратное ускорение в рамках одного GPU |77| и количество графических ядер в одном ускорителе увеличивается, то такая ситуация может быть достигнута в ближайшие пару .нет. Стоит отметить, что такая проблема в принципе невозможна па ускорителях Intel Xeon Phi.
Использование равномерной сетки в декартовых координатах дня решения уравнений гидродинамики позволяет использовать произвольную декартову топологию для декомпозиции расчетной области. Такая организация вычислений имеет потенциально бесконечную масштабируемость. В коде AstroPhi 2.0 используется многоуровневая одномерная декомпозиция расчетной области. По одной координате внешнее одномерное разрезание происходит средствами технологии MPI, внутри каждой подобласти разрезание происходит средствами ОрепМР, адаптированного дня MIC-архитектур (рис. 2). Такой подход использовался также и в первой версии программного кода AstroPhi |78| с учетом использования offload режима.
Такая декомпозиция связана с топологией и архитектурой гибридного суперЭВМ RSC PetaStream, который был использован для вычислительных экспериментов. Основные его параметры указаны в табл. 4.
Далее приведем исследования производительности различных подсистем кода.
3.1. Исследование ускорения кода в рамках одного MIC. Проводилось исследование ускорения кода AstroPhi на сетке 5123, для этого замерялось время каждого этапа (Euler, Lagrange) численного метода в секундах, а затем вычислялась их сумма (Total) при различном числе используемых логических ядер (Threads). Ускорение Р (Speed-Up) вычислялось по формуле 1:
Total1
TotalK '
Таблица 4
Основные параметры суперкомпьютера RSC PetaStream
Характеристика Архитектура Вычислительные узлы Логическое количество узлов Логическое количество ядер Оперативная память Хранение данных
Значение MIC (режим native) Intel Xeon Phi 7120 D 64 15372
1008 ГБ/узел GDDR5 Intel SSD DC S3500/S3700
где Tota/i — время вычислений та одном логическом ядре, Tota/K — время вычислений при использовании K логических ядер. Также была сделана оценка реальной производительности (SVeaO- Для оценки пиковой производительноети Speafc ускорителя Intel Xeon Phi была использована формула 2:
где Units = 1FLOPS — число скалярных устройств, Frequency = 1,238GHz — тактовая частота, Cores = 60 — число физических ядер. Для нахождения эффективности использования ускорителя Smtio будем использовать формулу 3
Эта величина производительности с одной стороны несколько странная, так как пиковый уровень производительности Intel Xeon Phi — это величина порядка одного терафлоп-са. Однако — такая производительность в основном достигается векторным расширением, которые мы пока не используем. Результаты исследований ускорения и производительности на сетке 5123 приведены на рис, 3,
Таким образом, было получено 134-кратное ускорение (масштабируемость в сильном смысле) в рамках одного ускорителя Intel Xeon Phi, Количество гигафлопсов соотносится с результатами, полученными в книге [95], А результат порядка чуть менее 29 GFLOPS — это примерно 40 % от пиковой, что достаточно высокий результат,
3,2, Исследование масштабируемости кода. Проводилось исследование масштабируемости кода AstroPhi на расчетной сетке 512p х 512 х 512 при использовании 4 логических ядер на каждый ускоритель, где p — число используемых ускорителей. Таким образом, па каждый ускоритель приходится размер подобласти 5123, Для исследования масштабируемости замерялось время каждого этапа (Euler, Lagrange) численного метода в секундах, а затем вычислялась их сумма (Total) при различном числе используемых ускорителей
T
Speak = Units x Frequency x Cores = 74,28GFLOPS,
(2)
(3)
T
Tota/1
(4)
Tota/P '
'p
где То^а^ — время вычислений на одном ускорителе при использовании одного ускорителя, То£а/р — время вычислений на одном ускорителе при использовании р ускорителей. Результаты исследований ускорения приведены на рис, 4,
140 1201003 80 "О
CÛ 60 0 о.
со
40
200-
30-, 2520-W 1RJ
о. 15-1
0
1 ЮН О
5
0-1
2 4 8 16 32 64 128 256
The logical cores
(b)
2 4 8 16 32 64 128 256
The logical cores
(a)
Рис.3. Ускоренно (а) и реальная производительность (b) кода AstroPhi 2.0 в рамках одного Intel Xoori Phi (а)
1,2 1,0 0,8
3 0,6
го
го
я °'4
0,2-
0,0-ттт-1
16
" I 1
32
ггп-64
MICs
Рис. 4. Масштабируемость кода AstroPhi 2.0 при различном числе ускорителей Intel Xeon Phi
Таким образом, при использовании вссх 64 ускорителей была получена 92 % эффективность, что является достаточно высоким результатом.
3.3. Исследование скорости чтения/записи в оперативной памяти. Как уже говорилось ранее, основной операцией в прочих этапах численного метода является фактически чтение элемента из одного массива и после простой операции запись в новый массив. Дня исследования этих процедур мы сформулировали очень простую дня анализа задачу — задачу копирования в рамках одного ускорителя Intel Xeon Phi массива размером 512 мегабайт, что соответствует памяти, используемой дня размещения одного гидродинамического параметра на сетке 5123, После достаточно простой оптимизации кода была получена скорость копирования данных 1921 МБ/с.
3.4. Исследование скорости записи на диск. Важной операцией, не связанной напрямую с процессом вычислений, является операция записи па диск промежуточных результатов. Дня этого мы сохраняли па диск трехмерный массив достаточно большой размерности, при этом запись каждого элемента проводилась в текстовом режиме с помощью следующего кода:
Phi!
Phi,
ti
Рис. 5. Схема передачи данных между двумя ускорителями Intel Xeon Phi
t
2
t
3
t
4
Таблица 5
Время копирования памяти в рамках одного ускорителя
Конфигурация Время (с)
ti 0,6251
t2 0,4180
t3 0,2812
t4 0,2907
for (i 0: i XX: i--)
for (k = 0;k<:XY;k--)
for (1=0:1<XZ: 1--)
{
fprintf (font ,"%lf %lf %lf %lf\n",i*h,k*h,l*h,a|index(i,k,l)|);
}
В результате скорость записи результатов па диск составила 3,26 МБ/с.
3.5. Исследование времени передали данных между ускорителями. Для исследования скорости передачи данных исследовалась передача данных между двумя ускорителями Intel Xeon Phi. Осуществлялась пересылка 64 массивов по 8 мегабайт каждый. Такой размер массива соответствует одному слою перекрытия для расчетной сетки 10243, Схему передачи данных можно проиллюстрировать рис. 5.
Сведем в табл. 5 времена каждого этапа.
В результате время обмена 64 массивов размером 8 мегабайт составило 1,615 с. Таким образом, скорость передачи данных составляет 635 МБ/с, а общий процент во времени счета составляет несколько процентов, что соответствует полученной масштабируемости.
Рис. 6. Результаты моделирования на момент времени £ =120 миллионов лет. Столбцовая плотность в М©ре-2 газовой (а) и бесстолкновительной (Ь) компонент, молекулярного водорода (с). Средняя скорость процесса звездообразования в М©ре-2тут-1 — (с1)
4. Результаты вычислительных экспериментов. Будем моделировать столкновение двух галактик с массой M = 1013Мо, летящих со скоростью vcr = 800 км/с, каждая из которых задана двумя самогравитирующими сферическими облаками дня описания газовой и бесстолкновительной компонент с равновесным начальным распределением плотности, давления/тензора дисперсии скоростей. Галактики вращаются в противоположные стороны с дифференциальным вращением:
д Ф
Результаты моделирования представлены па рис. 6.
На рис. 6 представлены результаты распределения столбцевой плотности газовой и звездной компоненты. После процесса столкновения за фронтом ударной волны происходит активный рост скорости звездообразования, а также па месте будущей галактики образуется молекулярный водород.
Заключение, В статье представлен новый гидродинамический программный код AstroPhi 2,0 для численного моделирования астрофизических процессов на гибридных суперЭВМ, оснащенных ускорителями Intel Xeon Phi, Описаны детали параллельной реализации кода и элементы со-дизайна численного алгоритма, которые позволили сделать эффективную программную реализацию.
Проведены вычислительные эксперименты по исследованию производительности различных подсистем кода на кластере RSC PetaStream, Для кода AstroPhi 2,0 было получено 134-кратное ускорение (масштабируемость в сильном смысле) в рамках одного ускорителя Intel Xeon Phi, 92-процентная эффективность (масштабируемость в слабом смысле) при использовании 64 ускорителей Intel Xeon Phi, и 40-процентная производительность от пиковой при использовании скалярных вычислений на ускорителе Intel Xeon Phi,
Приведены результаты суперкомпьютерных вычислительных экспериментов по столкновению дисковых галактик. Смоделированы области с активным звездообразованием и образованием молекулярного водорода, которые получаются в результате столкновения галактик.
Авторы выражают благодарность коллегам из организаций Intel (Николаю Меетеру и Дмитрию Петунину) и RSC Group (Александру Московскому, Павлу Лавренко и Борису Гагаринову) за предоставление доступа к кластеру RSC PetaStream и подробные консультации по его использованию.
Список литературы
1. Kim J., Park С., Richard Gott III J., Dubinski J. The Horizon Run N-Bodv Simulation: Barvon Acoustic Oscillations and Topology of Large-scale Structure of the Universe // The Astrophvsical Journal. 2009. V. 701, I. 2. P. 1547-1559.
2. Springel V., et al. Simulations of the formation, evolution and clustering of galaxies and quasars // Nature. 2005. V. 435. P. 629-636.
3. Teyssier R., et al. Full-skv weak-lensing simulation with 70 billion particles // Astronomy k, Astrophysics. 2009. V. 497, I. 2. P. 335-341.
4. Klypin A., Trujillo-Gomez S., Primack J. Dark Matter Halos in the Standard Cosmological Model: Results from the Bolshoi Simulation // The Astrophvsical Journal. 2011. V. 740, I. 2, 102. P. 1 17.
5. Genel S., et al. Introducing the Illustris project: the evolution of galaxy populations across cosmic time // Monthly Notices of the Royal Astronomical Society. 2014. V. 445, I. 1. P. 175-200.
6. Vogelsberger XL. et al. Introducing the Illustris Project: simulating the coevolution of dark and visible matter in the Universe // Monthly Notices of the Royal Astronomical Society. 2014. V. 444, I. 2. P. 1518-1547.
7. Vogelsberger XL. et al. Properties of galaxies reproduced by a hydrodvnamic simulation // Nature. 2014. V. 509. P. 177-182.
8. Chain R., et al. The EAGLE simulations of galaxy formation: calibration of subgrid physics and model variations // Monthly Notices of the Royal Astronomical Society. 2015. V. 450, I. 2. P. 1937-1961.
9. Schaye J., et al. The EAGLE project: simulating the evolution and assembly of galaxies and their environments // Monthly Notices of the Royal Astronomical Society. 2015. V. 446, I. 1. P. 521554.
10. mllgrom m. A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis // The Astrophvsical Journal. 1983. V. 270. P. 365-370.
11. Chilingarian I., Di Matteo P., Combes F., Melchior A., Semelin B. The GalMer database: galaxy mergers in the virtual observatory // Astronomy k Astrophysics. 2010. V. 518, A61. P. 1-14.
12. 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. V. 55, I. 9. P. 770-783.
13. Schweizer F. Merger-Induced Starbursts // Astrophysics and Space Science Library. 2005. V. 329. P. 143-152.
14. Sol Alonso M., Lamb as D., Tissera P., Cold well G. Active galactic nuclei and galaxy interactions // Monthly Notices of the Royal Astronomical Society. 2007. V. 375, I. 3. P. 1017-1024.
15. Blecha L., Loeb A., Narayan R. Double-peaked narrow-line signatures of dual supermassive black holes in galaxy merger simulations // Monthly Notices of the Royal Astronomical Society. 2013. V. 429, I. 3. P. 2594-2616.
16. Rodriguez C., Taylor G., Zavala R., Pihlstrom Y., Peck A. Hi observations of the supermassive binary black hole system in 0402+379 // The Astrophvsical Journal. 2009. V. 697, I. 1. P. 37-44.
17. Combes F., Melchior A. Chemodvnamical evolution of interacting galaxies // Astrophysics and Space Science. 2002. V. 281, I. 1-2. P. 383-387.
18. Glover S., Mac Low M. Simulating the Formation of Molecular Clouds. I. Slow Formation by Gravitational Collapse from Static Initial Conditions // The Astrophvsical Journal Supplement Series. 2007. V. 169, I. 2. P. 239-268.
19. Glover S., Mac Low M. Simulating the Formation of Molecular Clouds. II. Rapid Formation from Turbulent Initial Conditions // The Astrophvsical Journal. 2007. V. 659, I. 2. P. 1317-1337.
20. Ploeckinger S., Hensler G., Recchi S., Mitchell N., Kroupa P. Chemodvnamical evolution of tidal dwarf galaxies. I. Method and IMF dependence // Monthly Notices of the Royal Astronomical Society. 2014. V. 437, I. 4. P. 3980-3993.
21. Ploeckinger S., Recchi S., Hensler G., Kroupa P. Chemodvnamical evolution of tidal dwarf galaxies. II. The long-term evolution and influence of a tidal field // Monthly Notices of the Royal Astronomical Society. 2015. V. 447, I. 3. P. 2512-2525.
22. Recchi S. Chemodvnamical Simulations of Dwarf Galaxy Evolution // Advances in Astronomy. 2014. V. 2014, 750754. P. 1-30.
23. Few C., Courty S., Gibson B., Michel-Dansac L., Calura F. Chemodvnamics of a simulated disc galaxy: initial mass functions and Type la supernova progenitors // Monthly Notices of the Royal Astronomical Society. 2014. V. 444, I. 4. P. 3845-3862.
24. mlnchev I., Chiappini C., Martig M. Chemodvnamical evolution of the Milky Way disk.
I. The solar vicinity // Astronomy k Astrophysics. 2013. V. 558, A9. P. 1-25.
25. Minchev I., Chiappini C., Martig M. Chemodvnamical evolution of the Milky Way disk.
II. Variations with Galactic radius and height above the disk plane // Astronomy k Astrophysics. 2014. V. 572, A92. P. 1-19.
26. pllklngton K., et al. Metallicity gradients in disks. Do galaxies form inside-out? // Astronomy k Astrophysics. 2012. V. 540, A56. p. 1-12.
27. Recchi S., Spitoni E., Matteucci F., Lanfranchi G. The effect of differential galactic winds on the chemical evolution of galaxies // Astronomy k Astrophysics. 2008. V. 489, I. 2. P. 555-565.
28. Brusadin G., Matteucci F., Romano D. Modeling the chemical evolution of the Galaxy halo // Astronomy k Astrophysics. 2013. V. 554, A135. P. 1-10.
29. Harvey-Smith L., Cohen R. Discovery of large-scale methanol and hvdroxvl maser filaments in W3(OH) // Monthly Notices of the Royal Astronomical Society. 2006. V. 371, 1. 4. P. 1550-1558.
30. Glinskiy В., Kulikov I., Snytnikov A., Romanenko A., Chernykh I., Vshivkov V. Co-design of Parallel Numerical Methods for Plasma Physics and Astrophysics // Supercomputing frontiers and innovations. 2014. V. 1, I. 3. P. 88-98.
31. Mitchell N., Vorobyov E., Hensler G. Collisionless Stellar Hydrodynamics as an Efficient Alternative to N-bodv Methods // Monthly Notices of the Royal Astronomical Society. 2013. V. 428. P. 2674-2687.
32. glncold R. A., Monaghan J. J. Smoothed particle hydrodynamics - Theory and application to non-spherical stars // Monthly Notices of the Royal Astronomical Society. 1977. V. 181. P. 375-389.
33. Lucy L. B. A numerical approach to the testing of the fission hypothesis // The Astrophvsical Journal. 1977. V. 82. P. 1013-1024.
34. O'Shea В., Bryan G., Bordner J., Norman XL. Abel Т., Harkness 1С. Kritsuk A. Adaptive Mesh Refinement - Theory and Applications // Lectures Notes of Computer Science Engineering. 2005. V. 41. P. 341-350.
35. Hockney R. W., Eastwood J. W. Computer Simulation Using Particles. N.Y.: McGraw-Hill, 1981.
36. Couchman H.M.P. Mesh-refined P3M: A fast adaptive N-bodv algorithm // The Astrophvsical Journal. 1991. V. 368. P. L23-L26.
37. Barnes J., Hut P. A hierarchical O(NLogN) force-calculation algorithm // Nature. 1986. V. 324. P. 446-449.
38. Dubinski J., Kim J., Park C., Humble R. GOTPM: a parallel hybrid particle-mesh treecode // New Astronomy. 2004. V. 9. P. 111-126.
39. Fedorenko R. A relaxation method for solving elliptic difference equations // USSR Computational Mathematics k, Mathematical Physics. 1961. V. 1. P. 1092-1096.
40. Godunov S.K. A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodvnamic Equations // Mathematicheskv Sbornik (in russian). 1959. V. 47. P. 271-306.
41. Kulikovskii A.G., Pogorelov N.V., Semenov A.Yu. Mathematical Aspects of Numerical Solution of Hyperbolic Systems. XL: Fizmatlit (in russian), 2001.
42. toro E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Heidelberg: SpringerVerlag, 1999.
43. Courant R., Isaacson E., Rees M. On the solution of nonlinear hyperbolic differential equations by finite differences // Communications on Pure and Applied Mathematics. 1952. V. 5. P. 243-256.
44. Roe P. Approximate Riemann solvers, parameter vectors, and difference solvers // Journal of Computational Physics. 1997. V. 135. P. 250-258.
45. engquist В., Osher S.J. One-sided difference approximations for nonlinear conservation laws // Mathematics of Computational. 1981. V. 36. P. 321-351.
46. Harten A., Lax P.D., Van Leer B. On upstream differencing and Godunov-tvpe schemes for hyperbolic conservation laws // Society for Industrial and Applied Mathematics. 1983. V. 25. P. 35-61.
47. Einpeld B. On Godunov-tvpe methods for gas dynamics // SIAM Journal of Numerical Analysis. 1988. V. 25. P. 294-318.
48. Batten P., Clarke N., Lambert C., Causon D. M. On the Choice of Wavespeeds for the HLLC Riemann Solver // SIAM Journal of Computing. 1997. V. 18. P. 1553-1570.
49. Van Leer B. Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov's Method // Journal of Computational Physics. 1979. V. 32. P. 101-136.
50. Jin S., Xin Z. The Relaxation Schemes for Systems of Conservation Laws in Arbitrary Space Dimensions // Communications on Pure and Applied Mathematics. 1995. V. 48. P. 235-276.
51. Collela P., Woodward P.R. The Piecewise Parabolic Method (PPM) Gas-Dvnamical simulations // Journal of Computational Physics. 1984. V. 54. P. 174-201.
KyjiuKoe M. M., lIepHux M. R.
23
52. Liu X., Osher S., Chan T. Weighted essentially non-oscillatory schemes // Journal of Computational Physics. 1994. V. 115. P. 200-212.
53. Godunov S. K., Manuzina Yu. D., Nazareva \ LA. Experimental analysis of convergence of the numerical solution to a generalized solution in fluid dynamics // Computational Mathematics and Computational Physics. 2011. V. 51. P. 88-95.
54. Popov XL. Ustyugov S. Piecewise parabolic method on local stencil for gasdvnamic simulations // Computational Mathematics and Mathematical Physics. 2007. V. 47, I. 12. P. 19701989.
55. Popov XL. Ustyugov S. Piecewise parabolic method on a local stencil for ideal magnetohydrodvnamics // Computational Mathematics and Mathematical Physics. 2008. V. 48, I. 3. P. 477-499.
56. Pearcea F. R., Couchman H. M. P. Hydra: a parallel adaptive grid code // New Astronomy. 1997. V. 2. P. 411-427.
57. Wadsley J.W., Stadel J., Quinn T. Gasoline: a flexible, parallel implementation of TreeSPH // New Astronomy. 2004. V. 9. P. 137-158.
58. Matthias S. GRAPESPH: cosmological smoothed particle hydrodynamics simulations with the special-purpose hardware GRAPE // Monthly Notices of the Royal Astronomical Society. 1996. V. 278. P. 1005-1017.
59. Springel V. The cosmological simulation code GADGET-2 // Monthly Notices of the Royal Astronomical Society. 2005. V. 364. P. 1105-1134.
60. zlegler U. Self-gravitational adaptive mesh magnetohydrodvnamics with the NIRVANA code // Astronomy k Astrophysics. 2005. V. 435. P. 385-395.
61. Mignone A., Plewa T., Bodo G. The Piecewise Parabolic Method for Multidimensional Relativistic Fluid Dynamics // The Astrophvsical Journal. 2005. V. 160. P. 199-219.
62. Hayes J., Norman XL. Fiedler R. et al. Simulating Radiating and Magnetized Flows in Multiple Dimensions with ZEUS-MP // The Astrophvsical Journal Supplement Series. 2006. V. 165. P. 188-228.
63. Teyssier R. Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES // Astronomy k Astrophysics. 2002. V. 385. P. 337-364.
64. Kravtsov A., Klypin A., Hoffman Y. Constrained Simulations of the Real Universe. II. Observational Signatures of Intergalactic Gas in the Local Supercluster Region // The Astrophvsical Journal. 2002. V. 571. P. 563-575.
65. Stone J. et al. Athena: A New Code for Astrophvsical MHD // The Astrophvsical Journal Supplement Series. 2008. V. 178. P. 137-177.
66. Brandenburg A., Dobler W. Hvdromagnetic turbulence in computer simulations // Computer Physics Communications. 2002. V. 147. P. 471-475.
67. Gonzalez XL. Audit E., Huynh P. HERACLES: a three-dimensional radiation hydrodynamics code // Astronomy k Astrophysics. 2007. V. 464. P. 429-435.
68. Krumholz M.R., Klein R.I., McKee C.F., Bolstad, J. Equations and Algorithms for Mixed-frame Flux-limited Diffusion Radiation Hydrodynamics // The Astrophvsical Journal. 2007. V. 667. P. 626-643.
69. Mignone A. et al. PLUTO: a Numerical Code for Computational Astrophysics // The Astrophvsical Journal Supplement Series. 2007. V. 170. P. 228-242.
70. Almgren A. et al. CASTRO: A New Compressible Astrophvsical Solver. I. Hydrodynamics and Self-gravity // The Astrophvsical Journal. 2010. V. 715. P. 1221-1238.
71. Schive H., Tsai Y., Chiueh T. GAMER: a GPU-accelerated Adaptive-Mesh-Refinement Code for Astrophysics // The Astrophvsical Journal. 2010. V. 186. P. 457-484.
72. Murphy J., Burrows A. BETHE-Hvdro: An Arbitrary Lagrangian-Eulerian Multidimensional Hydrodynamics Code for Astrophvsical Simulations // The Astrophvsical Journal Supplement Series' 2008. V. 179. P. 209-241.
73. Springel V. E pur si muove: Galilean-invariant cosmological hydrodvnamical simulations on a moving mesh // Monthly Notices of the Royal Astronomical Society. 2010. V. 401. P. 791-851.
74. Bruenn S. et al. 2D and 3D core-collapse supernovae simulation results obtained with the CHIMERA code // Journal of Physics. 2009. V. 180. P. 1-5.
75. Hopkins P. A new class of accurate, mesh-free hydrodvnamic simulation methods // Monthly Notices of the Royal Astronomical Society. 2015. V. 450, I. 1. P. 53-110.
76. Vshivkov V., Lazareva G., Snytnikov A., Kulikov I., Tutukov A. Hydrodvnamical code for numerical simulation of the gas components of colliding galaxies // The Astrophvsical Journal Supplement Series. 2011. V. 194, I. 47. P. 1-12.
77. Kulikov I. GPUPEGAS: A New GPU-accelerated Hydrodvnamic Code for Numerical Simulations of Interacting Galaxies // The Astrophvsical Journal Supplements Series. 2014. V. 214, I. 12. P. 1-12.
78. Kulikov I.M., Chernykh I.G., Snytnikov A.V., Glinskiy B.M., Tutukov A.V. AstroPhi: A code for complex simulation of dynamics of astrophvsical objects using hybrid supercomputers // Computer Physics Communications. 2015. V. 186. P. 71-80.
79. Bergin E., Hartmann L., Raymond J., Ballesteros-Paredes J. Molecular Cloud Formation behind Shock Waves // The Astrophvsical Journal. 2004. V. 612. P. 921-939.
80. tlelens A., Hollenbach D. Photodissociation regions. I — Basic model. II — A model for the Orion photodissociation region // The Astrophvsical Journal. 1985. V. 291. P. 722-754.
81. Cazaux S., Tielens A. H2 Formation on Grain Surfaces // The Astrophvsical Journal. 2004. V. 604. P. 222-237.
82. Caselli P., Walmsley C., Terzieva R., Herbst E. The Ionization Fraction in Dense Cloud Cores // The Astrophvsical Journal. 1998. V. 499. P. 234-249.
83. Bergin E., Plume R., Williams J., Myers P. The Ionization Fraction in Dense Molecular Gas. II. Massive Cores // The Astrophvsical Journal. 1999. V. 512. P. 724-739.
84. Van Der Tak F., Van Dishoeck E. Limits on the cosmic-rav ionization rate toward massive young stars // Astronomy k, Astrophysics. 2000. V. 358. P. L79-L82.
85. Draine B., Bertoldi F. Structure of Stationary Photodissociation Fronts // The Astrophvsical Journal. 1996. V. 468. P. 269-289.
86. Draine B. Photoelectric heating of interstellar gas // The Astrophvsical Journal Supplement Series. 1978. V. 36. P. 595-619.
87. Katz N., Weinberg D., Hernquist L. Cosmological simulations with TreeSPH // The Astrophvsical Journal Supplement Series. 1996. V. 105. P. 19-35.
88. Springel V., Hernquist L. Cosmological smoothed particle hydrodynamics simulations: a hybrid multiphase model for star formation // Monthly Notices of the Royal Astronomical Society. 2003. V. 339, I. 2. P. 289-311.
89. Sutherland R., Dopita M. Cooling functions for low-density astrophvsical plasmas // The Astrophvsical Journal Supplement Series. 1993. V. 88. P. 253-327.
90. Kulikov I., Vorobyov E. Using the PPML approach for constructing a low-dissipation, operator-splitting scheme for numerical simulations of hydrodvnamic flows // New Astronomy. 2015. (in progress).
91. Vshivkov V., Lazareva G., Snytnikov A., Kulikov I., Tutukov A. Computational methods for ill-posed problems of gravitational gasodvnamics // Journal of Inverse and Ill-posed Problems. 2011. V. 19, I. 1. P. 151-166.
92. Godunov S., Kulikov I. Computation of Discontinuous Solutions of Fluid Dynamics Equations with Entropy Nondecrease Guarantee // Computational Mathematics and Mathematical Physics. 2014. V. 54, I. 6. P. 1012-1024.
93. Frigo M., Johnson S. The Design and Implementation of FFTW3 // Proceedings of the IEEE. 2005. V. 93, I. 2. P. 216-231.
94. Kalinkin A., Laevsky Y., Gololobov S. 2D Fast Poisson Solver for High-Performance Computing // Lecture Notes in Computer Science. 2009. V. 5698. P. 112-120.
95. Jeppers J., Reinders J. Intel Xeon Phi Coprocessor High Performance Programming. Elsevier, 2013.
Куликов Игорь Михайлович — канд. физ.-мат. наук, науч. сотр. Института вычислительной математики и математической геофизики СО РАН, e-mail: kulikov@ssd.sscc.ru.
Черных Игорь Геннадьевич — канд. физ.-мат. наук, старш. науч. сотр. Института вычислительной математики и математической геофизики СО РАН, e-mail: chernykh@parbz.sscc.ru.
Дата поступления — 28.06.2015
