CC BY
50
Известия Томского политехнического университета. Инжиниринг георесурсов. 2015. Т. 326. № 8
УДК 519.6
PADME - НОВЫЙ КОД ДЛЯ МОДЕЛИРОВАНИЯ ПРОЦЕССА ФОРМИРОВАНИЯ ГЕОРЕСУРСОВ ПЛАНЕТ НА ГЕТЕРОГЕННЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ
Протасов Виктор Александрович,
магистрант факультета Прикладной математики и информатики Новосибирского государственного технического университета, Россия, 630073, г. Новосибирск, пр. Карла Маркса, 20. E-mail: [email protected]
Куликов Игорь Михайлович,
канд. физ.-мат. наук, мл. науч. сотр. лаборатории Параллельных алгоритмов решения больших задач Института Вычислительной математики и математической геофизики Сибирского отделения Российской академии наук, Россия, 630090, г. Новосибирск, пр. Академика Лаврентьева, 6. E-mail: [email protected]
Актуальность работы. В последние годы было обнаружено большое количество планет, но возникают сложности с объяснением того, как они образуются. До недавнего времени единственным объектом для наблюдения являлась Солнечная система, и все гипотезы формирования планет основывались именно на этих наблюдениях. Со временем образовалось достаточно четкое понимание того, как происходило формирование Солнечной системы, но, тем не менее, остаются некоторые сомнения, так как неизвестно, что было в начале этого процесса, а что было приобретено позднее. Кроме того, сформировавшиеся представления зачастую не могут объяснить особенностей других планетных систем. Также очень важен вопрос поиска планет земного типа. Даже если какая-то экзопланета будет обладать схожими с Землей характеристиками, нельзя однозначно утверждать, что мы нашли «вторую Землю», так как внутренний, геологический, состав может существенно отличаться. Яркий пример тому - Венера. Другой актуальный вопрос - освоение георесурсов не только соседних планет и астероидов, но в скором времени и более удаленных. Если заранее удастся узнать какими полезными ископаемыми обладает тот или иной космический объект, то можно существенно сократить расходы на анализ грунта с помощью дорогостоящих космических аппаратов. Поэтому важно знать не только кинематические параметры планет, но и их внутренний состав. Таким образом, актуальным становится вопрос о моделировании химокинетических процессов уже на ранних стадиях эволюции планетной системы. Математическое моделирование наравне с наблюдением может помочь найти ответы на этот вопрос. Но вычислительная астрофизика, как и многие другие области науки, очень требовательна к ресурсам компьютерных систем, если необходимо получить высококачественное решение. Поэтому вопрос разработки новых численных методов и математических моделей также актуален, как и более эффективное использование имеющихся вычислительных мощностей для уже существующих методов.
Цель исследования: разработка нового метода для моделирования процесса планетообразования в 3D2V постановке на основе двухфазного подхода, адаптированного для использования в гетерогенных вычислительных системах, оснащенных графическими ускорителями с поддержкой технологии NVIDIA CUDA.
Методы исследования. Для моделирования газовой компоненты используется метод крупных частиц Белоцерковского-Давыдова, модифицированный с использованием метода Годунова. Пылевая компонента описывается системой N тел, динамика которой просчитывается Particle-Mesh методом. Для повышения точности моделирования динамики частиц используется подход Clouds-in-Cells. Уравнение Пуассона для гравитационного потенциала решается методом быстрого преобразования Фурье. Результаты. Разработан новый метод для моделирования процесса планетообразования. Представлены результаты тестирования. Газодинамическая часть была проверена на модельных задачах газовой динамики, а оценка правильности решения уравнения Пуассона выполнена на функции с известным распределением потенциала. Также приведен результат моделирования газопылевого диска с образованием уплотнения из газа и пыли, которое можно интерпретировать как протопланету. Показана целесообразность использования графических ускорителей для такого рода задач.
Ключевые слова:
Математическое моделирование, вычислительная астрофизика, гравитационная газовая динамика, планетообразование, внутреннее строение планет, параллельные численные методы, гетерогенные вычислительные системы, GPGPU, CUDA.
Введение
Поиск планетных систем - одна из важнейших задач современной астрономии. Более десятилетия назад единственным объектом для изучения была Солнечная система. Факт существования других планетных систем был недоказуем, ввиду отсутствия достаточных для этого возможностей. Но в начале 90-х гг., с появлением достаточно мощных телескопов, были обнаружены первые внегалактические планеты, число которых к концу 2014 г. стало порядка 1800.
Современные представления о возникновении Солнечной системы основаны на гипотезе Канта-Лапласа, в соответствии с которой Солнце формировалось одновременно с планетами, которые образовались из околозвездного газа и пыли. Горячие компоненты во внутренней части диска образовали планеты земной группы, а сгустки вещества в холодной внешней части диска образовали ядра планет-гигантов, которые притягивали к себе частицы пыли и остывшего газа, образовавших в дальнейшем оболочку планет. Эта теория доста-
61
Протасов В.А., Куликов И.М. PADME - новый код для моделирования процесса формирования георесурсов ... С. 61-70
точно точно описывает структуру Солнечной системы, но дает сбой на некоторых других системах. Например, были обнаружены горячие планеты-гиганты, вращающиеся слишком близко к звезде. Также остается открытым вопрос о времени, необходимом для формирования планетной системы. Предполагалось, что этот процесс занимает сотни миллионов лет, но на деле это значение оказалось в десятки раз меньше.
Другая не менее важная проблема - формирование планет в двухзвездных системах. Приблизительно 20 % обнаруженных планет находятся именно в таких системах, но до сих пор нет однозначного объяснения их появления. Наблюдательные возможности позволяют обнаружить планеты в системах со звездами массой M«M0. Звезды с меньшей массой слишком тусклые, чтобы можно было изучить их в деталях, а звезды с большей массой настолько яркие, что даже ближайшие планеты не создают какого-либо значительного затенения. Кроме того, высокие скорости вращения таких звезд исключают возможность использования спектральных методов поиска. Все это делает поиск планетных систем очень трудоемким. Изучение обнаруженных планетных систем позволило описать несколько сценариев их образования. Например, в системе из двух звезд, одна из которых красный гигант, а другая белый карлик, может происходить сброс вещества первой звездой, часть которого затем притягивается к белому карлику, образуя протопланетный диск. Кроме этого, планетные системы могут образоваться в ходе эволюции контактных двойных звезд. Более подробно эти и некоторые другие теории рассмотрены в [1].
Все большую популярность набирают исследования геофизических свойств экзопланет, так как ответ на этот вопрос может помочь не только в поиске планет, схожих с Землей, но и позволит лучше понять процессы, происходящие в недрах планет. Конечно, наилучшим методом исследования является эксперимент, но даже для самых простых случаев - водородных планет - невозможно воссоздать тот диапазон давлений и температур, которые наблюдаются в реальных ситуациях. По наблюдениям также довольно сложно дать ответ на поставленный вопрос, так как даже определение массы и радиуса планет не всегда удается выполнить с достаточной точностью. Но если это все же удалось сделать, можно с некоторой долей вероятности определить некоторые параметры ядра и мантии экзопланеты [2]. Несмотря на то, что физика поведения вещества в недрах планет очень сложна, моделирование этого процесса остается возможным.
Запасы ресурсов нашей планеты сильно ограничены, а за последние годы объемы потребления и добычи основных источников энергии только возрастают [3]. Поэтому вопрос разработки внеземных месторождений не менее актуален, чем развитие альтернативных источников энергии, хоть и рассчитан на более удаленное будущее. Знание то-
го, какими свойствами обладают космические объекты - планеты, их спутники, астероиды и т. д., содержащие в себе интересующие нас минералы, - позволит исключить заведомо бесперспективные исследования грунта с использованием дорогостоящей аппаратуры, которую к тому же вряд ли удастся использовать повторно.
На сегодняшний день выделяют два основных подхода к моделированию протопланетного диска - двухкомпонентный подход [4] и двухфазный [5, 6]. Их отличие состоит в том, как моделируется пылевая компонента диска. В первом подходе частицы пыли представляются несжимаемой сплошной средой (жидкими частицами). Это позволяет учитывать обмен энергией и импульсом между газом и частицами. Во втором подходе частицы моделируются набором дискретных тел, которые оказывают влияние на гравитационное поле диска, но при этом не учитывается их взаимодействие с газом.
В данной работе представлен двухфазный подход. Для моделирования газовой компоненты диска применяется метод крупных частиц Белоцерковского-Давыдова [7], который является развитием метода частиц в ячейках [8], модифицированный с использованием схемы Годунова [9]. Частицы пыли представлены системой N-тел, для решения которой используется Particle-Mesh метод [10, 11]. Для вычисления гравитационного взаимодействия частиц и газа решается уравнение Пуассона для гравитационного потенциала с использованием FFT.
Постановка задачи
Рассмотрим модель динамики протопланетного диска в декартовых координатах, которая описывается следующей системой уравнений:
dp + div(pv) = 0, (1)
—- + div (vpv) = - grad (p) - pgrad (Ф), (2)
dt
+ div(psv) = -(y - l)psdiv(v), (3)
dt
dpE + div(pEv) = -div(pv) - (pgrad (Ф)^), (4)
dt
div( grad (Ф)) = 4nG p, (5)
P = (y- 1)pL, (6)
и pv2 pE = pL+—, (7)
d2 x F
dt2 m ’ (8)
F = -grad (Ф). (9)
62
Известия Томского политехнического университета. Инжиниринг георесурсов. 2015. Т. 326. № 8
Уравнения (1)-(7) описывают динамику самогра-витирующего газа, а уравнения (8), (9) - динамику произвольной частицы. Здесь р - плотность газа; v -скорость движения газа; p - давление газа; s - внутренняя энергия газа; E - полная энергия газа; у -показатель адиабаты; Ф - гравитационный потенциал; G - гравитационная постоянная; р - плотность распределения газа и частиц в трехмерной области; x - позиция некоторой частицы; m - ее масса.
Процесс планетообразования является существенно трехмерным. Но, с течением времени, изначально сферический объект под действием сил притяжения сжимается в центральной части и растягивается под действием центробежной силы. Поэтому, начиная с некоторого момента времени, он может быть рассмотрен в двумерном приближении. Но есть один нюанс. В трехмерном случае решение уравнения Пуассона для гравитационного
потенциала пропорционально —, тогда как для
r
двумерного случая - ln(—^ . Учитывая, что процесс все же является трехмерным, будем рассматривать уравнение Пуассона в трехмерной постановке для получения более физичного решения, а динамику газа и движения частиц - в двумерной.
Исходная задача неустойчива, а ее постановка некорректна по Адамару, поэтому необходим дополнительный механизм для контроля правильности получаемого решения [12]. Для этого будем проверять выполнение законов сохранения массы, импульса и энергии.
Моделирование газовой динамики
Моделирование газовой компоненты осуществляется с использованием метода крупных частиц Белоцерковского-Давыдова. Он является развитием метода Харлоу, который используется, например, в [13]. В его основе лежит схема расщепления по физическим процессам - решение исходной системы разбивается на два этапа:
1. Эйлеров этап, на котором решается система уравнений газовой динамики без адвективных членов, т. е. происходит пересчет параметров в предположении, что газ является неподвижным, с постоянной плотностью.
2. Лагранжев этап, на котором происходит адвективный перенос газодинамических величин.
В расчетной области вводится прямоугольная сетка с координатами узлов
х,. = х0 + ip. i = 0, Nx. yj = y0 + jhy . j = 0, Ny .
Все газодинамические параметры задаются в центрах ячеек:
Для конечно-разностной аппроксимации уравнений эйлерова этапа используется схема Годунова, которая для уравнений вида
dU (х) dF (x,U)
dt дх
записывается как
U”+l - U" At
F — - F —
i+— i —
----2------2 = 0,
Ax
где F — - значение на границе ячейки справа и
_ 2
слева. Для определения этих значений рассмотрим следующую систему уравнений газовой динамики для идеального газа на эйлеровом этапе:
dv 1 др
dt р дх’
др
~t
dv
(У- —) Р .
дх
Эта система имеет аналитическое решение:
v, = VL+VR + Pl - pR
Pl + Pr
PlPr (у- 1)(Pl + Pr )
• Pl+Pr + VL - VR Pl Pr (у- 1)( Pl + Pr )
2 2 V Pl + Pr ’
где (vL,pL,pL) - значения газодинамических параметров в левой ячейке, а (vR,pR,pR) - в правой. Таким образом, конечно-разностная схема решения системы уравнений на эйлеровом этапе имеет вид:
^и+1
Pjj -Pj At
= 0,
pvx - pvx
At
p - p ф" - Ф"
fi+1,j t'i-1,j n-1 i +1j i -1j
■P,„
2h
(10)
pv" -pv"
И y.,j И y.,j =
At
Кj+— - Pij-i - „-— % +i -ф"j -i
h PiJ 2h
(11)
= -(у-1) Р"
pn. - pn-1
At
( v* - vх
‘+—j___i-—.
h
v - v
у,., +i у,,
(12)
Уравнения (10)и (11) получены из уравнения (2) и позволяют вычислить х- и р-компонен-ту импульса, соответственно. Уравнение для давления (12) получено из (3), если в него подставить (6).
Для определения схемы расчета на лагранже-вом этапе рассмотрим произвольную ячейку (рис. 1). В этой ячейке в общем случае есть как исходящий поток газа, так и входящий. Будем рассматривать потоки только через ребра ячейки.
63
Протасов В.А., Куликов И.М. PADME - новый код для моделирования процесса формирования георесурсов ... С. 61-70
Рис. 1. Поток газа через ячейку Fig. 1. Gas flow through a cell
Количество вещества, которое останется в ячейке через время т, можно рассчитать по следующей формуле:
МП1 = Mf +-----(х _ flow _ inij - x _ flow _ out j ) +
x
+-----(y flow in.
h - ,,J
У _ flow _ouf . ).
Для определения потоков через границы используется модификация классического способа расчета, которая учитывает возможный «скос» границы (рис. 2) за счет различных скоростей в узлах.
Рис. 2. Скос границы ячейки Fig. 2. Downwash of a cell boundary
Выходящий поток по оси x тогда рассчитывается как
х _ flow _ outi.
1
2
t V:
k=0
- j+l
M,,j, Vx
> 0
Mt+y, Vx
< 0
Метод крупных частиц Белоцерковского-Давыдова, как и многие другие методы, использующие явные конечно-разностные схемы, является условно устойчивым. Условие устойчивости заключается в том, чтобы CFL<1. Где CFL - число Куран-та-Фредриксона-Леви, которое определяется как
ПТСТ T(Cs + max(l Vx 1,1 Vy l))
CrL =-------------------,
min(hx, hy)
где cs =
- скорость звука в веществе.
Моделирование пылевой компоненты
При моделировании динамики частиц основная трудность состоит в определении силы, действующей на некоторую частицу со стороны других частиц. Кроме этого, необходимо учитывать взаимодействие гравитационных полей газа и частиц.
В данной работе используется метод Particle-Mesh, который позволяет сократить вычислительные затраты на этом этапе моделирования. Суть метода заключается в том, что мы разбиваем расчетную область на конечное множество ячеек и в каждой из них рассчитываем плотность частиц. Затем для полученного распределения плотности решается уравнение Пуассона для гравитационного потенциала. Зная потенциал можно без особого труда посчитать силу притяжения как F=-grad (Ф), где Ф вычисляется из уравнения (5) с учетом
того, что р Pgas+Pparticks.
У данного метода существует весомый недостаток - точность. При использовании классического подхода (когда сила, действующая на частицу, вы-„ mm.
числяется как F, = tG —'~pL) точность расчета
i*j rij
силы притяжения зависит только от точности суммирования. В этом же методе существует несколько источников погрешности:
• вычисление плотности частицы в ячейке. Как и в методе Харлоу, здесь не избежать колебаний плотности при переходе частицы из ячейки в ячейку, что приводит к нефизичным флуктуациям решения;
• вычисление гравитационной силы. Гравитационный потенциал привязывается не к частице, а к центру ячейки.
Чтобы уменьшить влияние этих факторов плотность частицы в ячейке и сила, действующая на нее, вычисляются по методу Clouds-In-Cells (CIC) [14]. При таком подходе считается, что координаты частицы - координаты центра массы «облака» конечного размера. Плотность такого облака распределяется между ячейками, в которые оно попало (рис. 3).
Рис. 3. Плотность частицы распределяется между ячейками
(i,j), (i+1,j), (i,j+1), (i+1,j+1)
Fig. 3. Particle density is distributed among the cells (i,j), (i+1,j), (i,j+1), (i+1,j+1)
Таким образом, плотность частиц в некоторой ячейке и сила, действующая на частицу, вычисляются как:
Р>р = t Wi j(x> У)Pc(x У^
ceclouds
Fc =-tWi,j(x> y)grad (Ф) j,
i, i
где
64
Известия Томского политехнического университета. Инжиниринг георесурсов. 2015. Т. 326. № 8
Щ, ( х, у) =
( | х - х | ^
I 1 —-------- I х
^ к )
( I у - у,-О
х | 1 - j |
I К )
1 х - х < кх&| у - у, < К
0, eia^a
Данный подход, конечно же, не решает проблемы полностью, но он позволяет существенно сократить ошибку вычислений, что подробно рассмотрено в [14].
Решение уравнения Пуассона
Для определения гравитационного потенциала в расчетной области необходимо решить уравнение Пуассона, которое с учетом обезразмеривания выглядит как
ДФ = 4пр.
Для этого зададим трехмерную декартову сетку, центральный слой которой совпадает с сеткой для моделирования динамики газа и частиц, и будем аппроксимировать оператор Лапласа 27-точечным разностным шаблоном. Представим потенциал и плотность в виде суперпозиции по собственным функциям оператора, в результате чего получим следующую схему решения уравнения в пространстве гармоник:
V,
4пР,,
( 2 . 2 ( nj^
1 ч1 - э- I i) j
(
,( nm И
х | 1 - 2sin21 Пт | | х
l 3 l Ny))
(, 2 . 2 ( nn^
х I 1 — sin21 — | |
l 3 1N ))
(13)
Ф(Г)| =- m - S + ^ + Sz - 1х + 1у + 1 ~ 310 ^
G И И2 2 \r |3
где m - общая масса газа и частиц;
r = V(х - х0)2 + (у - у0)2 + (z - z0)2
- расстояние от центра масс до точки с координатами (x,y,z); Sx, Sy и Sz - статические моменты инерции; Ix, Iy и Iz - осевые моменты инерции, а I0 - момент инерции, проведенный через точку (x,y,z) и точку начала координат.
Тестирование
Для верификации программной реализации выполнялось тестирование отдельно газодинамической части с использованием тестов Годунова [16] и на модельных задачах. Отдельно был протестирован метод решения уравнения Пуассона на задачах с аналитическим решением. Выполнено тестовое моделирование динамики частиц.
Приведем некоторые результаты расчетов модельных задач газовой динамики.
1. Неустойчивость Рэлея-Тейлора. Рассмотрим область [-0,5;0,5]х[-0,5;0,5], заполненную газом, находящемся в поле силы тяжести, с плот-
ностью р =
1, у < 0
2, у > 0
и давлением p=2,5-py. На
2.
границе раздела сред задано возмущение:
v}=A(y)[1+cos(6nx)], А(у)
10-2, |у| < 0,01 0, \у\ > 0,01.
Результат представлен на рис. 5, а. Неустойчивость Кельвина-Гельмгольца. Неустойчивость возникает, когда контактирующие среды движутся с различными скоростями и на границах раздела имеются возмущения. Зададим в области [-0,5;0,5]х[-0,5;0,5] газ с плот-
12, |у| < 0,25
ностью р = ( , , , давлением p=2,5 и ско-
[1, |у| > 0,25 '
Таким образом, решение уравнение Пуассона можно разбить на три этапа:
1. Перевод в пространство гармоник 4npjk.
2. Вычисление в нем значений потенциала по формуле 13).
3. Преобразование найденных значений обратно в функциональное пространство
Для перехода в пространство гармоник и обратно используется быстрое преобразование Фурье.
Краевые условия уравнения Пуассона полностью определяют решение задачи, поэтому их постановка является достаточно важной проблемой. Известно, что на бесконечном удалении от объекта можно считать потенциал нулевым, но обстоятельства складываются таким образом, что у нас нет возможности удовлетворить необходимым условиям. В работе [15] используется способ вычисления значений на границе через моменты инерции:
ростью vx =
-0,5, |у| < 0,25 0,5, |у| > 0,25
v}-=A(y)[1+cos(8nx)],
0,01, |у - 0,25 < 0,1
где А(у) = <
0,01, |у + 0,25|<0,1.
0, иначе
Результат представлен на рис. 5, б.
3. Неустойчивость Рихтмайера-Мешкова. Данный вид неустойчивости возникает на границе раздела двух сред при прохождении через нее ударной волны. Рассмотрим конфигурацию, представленную на рис. 4. В области задано два слоя газа различной плотности и ударная волна. Результат моделирования показан на рис. 6, а. Похожую картину можно наблюдать и в случае, когда область разрежения находится внутри
65
Протасов В.А., Куликов И.М. PADME - новый код для моделирования процесса формирования георесурсов ... С. 61-70
более плотного слоя газа [17, 18]. Пусть область [-1;1]х[-0,5;0,5] заполнена покоящимся газом с параметрами p=p=1. В области [0,375;0,625]х х[-0,125;0,125] находится разреженный газ с плотностью p=1/29 и давлением p=1. При х>0,75 задана ударная волна с параметрами p=3/2 и p=4/3. Результат расчета представлен на рис. 6, б.
Рис. 4. Начальные параметры газа для неустойчивости Рихт-майера~Мешкова
Fig. 4. Initial parameters of gas for Richtmyer-Meshkov insta-
bility
Проверим правильность решения уравнения Пуассона на задаче с известной функцией потенциала на вложенных сетках:
Ф( г)
Ъп , 2п о
---г +— г
5 Ъ
4п
15г
г > 1
Р( г)
2г3 - 2г2 +1, г < 1 0, г > 1
Как видно из таблицы, имеет место 2-й порядок аппроксимации.
Таблица. Результаты решения уравнения Пуассона Table. Results of solution of Poisson equation
Размер сетки Grid size Относительная норма погрешности Relative norm of error
163 6,425461e-003
323 1,593152e-003
643 4,012327e-004
1283 1,040856e-004
2563 2,462044e-005
Рис. 5. Распределение плотности газа для неустойчивости a) Рэлея~Тейлора при t=6,2; б) Кельвина-Гельмгольца при t=3,0 Fig. 5. Gas density distribution for instability of: a) Rayleigh~Taylor at t=6,2; b) Kelvin~Helmholtz at t=3,0
Рис. 6. Распределение плотности газа для неустойчивости Рихтмайера-Мешкова при: а) t=3,0; б) t=0,7 Fig. 6. Gas density distribution for Richtmyer-Meshkov instability at: a) t=3,0; b) t=0,7
66
Известия Томского политехнического университета. Инжиниринг георесурсов. 2015. Т. 326. № 8
Рис. 7. Ускорение решения уравнения Пуассона Fig. 7. Acceleration of Poisson equation solution
Параллельная реализация
Распараллеливание выполняется с использованием графических ускорителей с поддержкой технологии NVIDIA CUDA [19, 20], впрочем, использование MIC-архитектур также является возможным [21]. Самым долгим этапом моделирования является решение уравнения Пуассона. Часть работы по распараллеливанию уже проведена за нас - преобразование Фурье реализовано на GPU в библиотеке cuFFT. Ускорить этап можно, если перенести на видеокарту учет краевых условий и вычисление потенциала в пространстве гармоник, что позволит дополнительно сократить количество пересылок данных. В результате удалось добиться почти 10-кратного ускорения на довольно слабой мобильной видеокарте (рис. 7).
Перенос газодинамики на видеокарту, как показала практика, не является целесообразным. Для обеспечения контроля точности решения проверяются законы сохранения энергии, массы и импульса, что требует редукции всех газодинамических параметров. Скорость выполнения редукции на GPU становится сравнимой с однопоточной версией на CPU только для очень больших объемов данных. Если оставить выполнение редукции на процессоре, а другую часть перенести на видеокарту, то значительно возрастает количество операций копирования данных с устройства, что составляет примерно половину времени расчета на одном ядре процессора. Поэтому было решено выполнять газодинамический расчет только на процессоре.
CPU GPU
J Poisson Stage
! GasDynamics ParticleModelling V time
Рис. 8. Схема выполнения этапов моделирования Fig. 8. Diagram of modeling stage implementation
Моделирование динамики частиц в текущей реализации прекрасно параллелится - было получено 4-кратное ускорение при моделировании динамики 106 частиц относительно одного ядра процессора. Кроме того, этот этап может выполняться независимо от этапа моделирования динамики газа, поэтому можно предложить следующую схему выполнения, представленную на рис. 8.
Моделирование динамики газопылевого диска
Рассмотрим область [-2,5;2,5]х[-2,5;2,5], в которой задан газ с параметрами:
12, r < 2 fl, r < 2
p(r) = 1 , p(r) = 1 ,
[3, r < 0,25 [2, r < 0,25
m(r) = 2(2 - r)2.
Также в области Pe[0,5;2,0] равномерно распределено 50 000 частиц с массой т=5-105и угловой скоростью o(r)=2(2-r)2+s, где s - небольшое отклонение.
На рис. 9, а, б видно образование кольца из газа и пыли, в котором формируется уплотнение, состоящее в основном из частиц. Этот сгусток плотности можно интерпретировать как потенциальную планету.
Заключение
В работе представлен новый метод моделирования процесса планетообразования на основе метода крупных частиц и PM-метода, ориентированный на использование в гибридных вычислительных системах. Продемонстрирована работоспособность алгоритма на модельных задачах. Проведен расчет динамики газопылевого диска.
Представленный метод не является финальной точкой. Это лишь робкий шаг на пути к высококачественным алгоритмам в области вычислительной астрофизики, позволяющим с высокой точностью дать ответ на важнейшие вопросы, связанные с изучением вселенной. Из возможных улучшений можно выделить добавление процесса коагуляции ча-
67
Протасов В.А., Куликов И.М. PADME - новый код для моделирования процесса формирования георесурсов ... С. 61-70
Рис. 9. Моделирование газопылевого диска: a) t=0,39; б) t=0,4 Fig. 9. Modeling of gas-and-dust disc: a) t=0,39; b) t=0,4
стиц, уменьшение влияния вычислительных эффектов при моделировании газовой динамики, использование менее требовательного ко времени метода решения уравнения Пуассона. Также актуальным остается использование суперкомпьютеров -рассмотренный метод имеет очень неплохие харак-
теристики в плане масштабируемости, как показано в [22].
Работа выполнена при финансовой поддержке РФФИ в рамках научных проектов № 15-01-00508 и №15-31-20150 молавед, а также поддержана грантом Президента РФ МК-6648.2015.9.
СПИСОК ЛИТЕРАТУРЫ
1. Tutukov A.V., Fedorova A.V. Formation of Planets during the Evolution of Single and Binary Stars // Astronomy Reports. -2012. - V. 54. - №4. - P. 305-314.
2. Can we Constrain Interior Structure of Rocky Exoplanets from Mass and Radius Measurements? / C. Dorn et al. // A&A. -
2015.- V. 577. - P. 12-31.
3. BP Statistical Review of World Energy. June 2015. 64th ed. URL: http://www.bp.com/content/dam/bp/pdf/Energy-econo-mics/statistical-review-2015/bp-statistical-review-of-world-en-ergy-2015-full-report.pdf (дата обращения: 03.03.2015)
4. A Two-Phase Code for Protoplanetary Disks / S. Inaba et al. // A&A. - 2005. - V. 431. - P. 365-379.
5. SEREN - a New SPH Code for Star and Planet Formation Simulations. Algorithms and Tests / D.A. Hubber et al. // A&A. -2011.- V. 529. URL: http://www.aanda.org/articles/aa/pdf/ 2011/05/aa14949-10.pdf (дата обращения: 03.03. 2015).
6. Численное решение трехмерных задач динамики самогравити-рующих многофазных систем / В.А. Вшивков и др. // Научный вестник НГТУ. - 2011. - № 3 (44). - С. 69-80.
7. Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике. - М.: Наука, 1982. - 392 с.
8. Harlow F.H. The Particle-in-Cell Method for Numerical Solution of Problems in Fluid Dynamics // Proceedings of Symposium in Applied Mathematics. - New York, 1963. - V. 15. - № 10. - P. 269-288.
9. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. - М.: Физматлит, 2001. - 608 с.
10. Hockney R.W., Eastwood J.W. Computer Simulation Using Particles. -New York: Adam Hilger, IOP Publication Ltd., 1989. - 540 с.
11. Moscardini L. Cosmological Particle-Mesh N-Body Simulations and the Role of Initial Conditions // Rend. Sem. Mat. Univ. Pol. Torino. - 1993 - V. 51. - № 3. - P. 249-266.
12. Computational Methods for Ill-Posed Problems of Gravitational Gasodynamics / Vshivkov V.A. et al. // Journal of Inverse and Ill-Posed Problems. - 2011. - V. 19. - № 1. - P. 151-166.
13. Боронина М.А., Вшивков В.А., Дудникова Г.И. Неявная схема для решения уравнений Максвелла в областях с различными масштабами // Доклады Академии наук высшей школы Российской Федерации. - 2014. - № 4. - С. 39-46.
14. Birdsall C.K., Fuss D. Clouds-in-Clouds, Clouds-in-Cells Physics for Many-Body Plasma Simulation // Journal of Computational Physics. - 1997. - V. 135. - P. 141-148.
15. Параллельная реализация на суперЭВМ модели газовой компоненты самогравитирующего протопланетного диска / Вшивков В.А. и др. // Вычислительные технологии. - 2007. -Т. 12. - № 3. - С. 38-52.
16. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. - Berlin: Heidelberg; New York: Springer-Verlag, 1999. - 640 p.
17. Рыбакин Б.П., Шидер Н.И. Построение параллельных алгоритмов для решения задач гравитационной газовой динамики // Вычислительные методы и программирование. - 2010. - Т. 11. - С. 388-394.
18. Рыбакин Б.П. Численное моделирование задач газовой динамики на многопроцессорных ЭВМ // Параллельные вычисления и задачи управления: Доклады пятой Междунар. конф. -Москва, 2010. - С. 161-170.
19. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA. - М.: ДМК Пресс, 2010. - 232 с.
20. Сандерс Дж., Кэндрот Э. Технология CUDA в примерах. Введение в программирование графических процессоров. - М.: ДМК Пресс, 2011. - 232 с.
21. Куликов И.М., Черных И.Г., Глинский Б.М. AstroPhi: программный комплекс для моделирования динамики астрофизических объектов на гибридных суперЭВМ, оснащенных ускорителями Intel Xeon Phi // Вестник Южно-Уральского государственного университета. Серия: Вычислительная математика и информатика. - 2013. - Т. 2. - № 4. - С. 57-79.
22. Параллельная реализация численной модели столкновения галактик / И.М. Куликов и др. // Вестник Новосибирского государственного университета. Серия: Информационные технологии. - 2011. - Т. 9. - №4. - С. 71-78.
68
Поступила 11.03.2015 г.
Известия Томского политехнического университета. Инжиниринг георесурсов. 2015. Т. 326. № 8
UDC 519.6
PADME - A NEW CODE FOR MODELING PLANET GEORESOURCES FORMATION ON HETEROGENEOUS COMPUTING SYSTEMS
Viktor A. Protasov,
Novosibirsk State Technical University, 20, Karl Marks avenue, Novosibirsk,
630073, Russia. E-mail: [email protected]
Igor M. Kulikov,
Institute of Computational Mathematics and Mathematical Geophysics SB RAS, 6, Lavrentyev avenue, Novosibirsk, 630090, Russia. E-mail: [email protected]
Relevance of the research. Many planets were detected in last few years, but there is no clear understanding of how they are formed. The Solar system was the only object for observation until recently, and all hypotheses about planet formation were based only on it. The fairly clear understanding about Solar system formation was founded with time, but there are some doubts yet, because we do not know what was in the beginning of the process, and what was acquired afterwards. Moreover, the formed ideas often could not explain some features of other systems. Searching for Earth-like terrestrial planets is another very important problem. Even if any of found exoplanets will be similar to Earth, we could not say that it is the «second Earth» exactly, because its internal, geological, composition could be different. Venus is a vivid example of this. Another relevant issue is exploring geo assets not only of nearby planets and asteroids, but in the near future of more distant space objects. If we know what minerals are inside them, we will cut the costs on the ground analysis with expensive spacecrafts. Therefore, it is very important to know not only the kinematic characteristics of a planet, but its internal composition as well. Thus, the issue of chemical-kinetics modeling in the early stages of the planetary system evolution becomes urgent. Mathematical modeling on a par with observation could help to find the answers to this question. But the computational astrophysics, as many other fields of science, is very demanding to resources of computing systems, if we want to obtain the high quality solution. So developing new numerical methods and mathematical models are as relevant as more efficient use of computational power in existing methods.
The aim of the study is to develop a new method for modeling planet formation in 3D2V formulation based on two-phase approach, adapted for using in heterogeneous computing systems equipped with graphics accelerators supporting NVIDIA CUDA technology.
The methods of the study. Fluids-in-cells method of Belotserkovskii~Davydov, modified with using the Godunov method, is used to model the gas component. A dust component is described by N-body system solved with Particle-Mesh method. The Clouds-in-Cells approach is used to increase the accuracy of modeling of the particles dynamics. Poisson equation for gravitational potential is solved with fast Fourier transform method.
The results. The authors have developed the method for modeling planet formation. The verification results are introduced. The gas dynamics part was tested with use of model problems in gas dynamics, and correctness of solution of Poisson equation was assessed by using function with known potential distribution. The paper introduces as well the gas-dust disk modeling results with formation of sealing of gas and dust, which can be interpreted as potential exoplanet. Advisability of using the graphics accelerators for such problems is demonstrated.
Key words:
Mathematical modeling, computational astrophysics, gravitational gas dynamics, planet formation, planetary interior, parallel numerical methods, heterogeneous computing systems, GPGPU, CUDA
The research was financially supported by RFBR within the scientific projects no. 15-01-00508 and 15-31-20150 мол а вед, as well as by the grant of the President of the RF МК-6648.2015.9.
REFERENCES
1. Tutukov A.V., Fedorova A.V. Formation of Planets during the Evolution of Single and Binary Stars. Astronomy Reports, 2012, vol. 54, no. 4, pp. 305-314.
2. Dorn C. Can we Constrain Interior Structure of Rocky Exoplanets from Mass and Radius Measurements? A&A, 2015, vol. 577, pp. 12-31.
3. BP Statistical Review of World Energy. June 2015. 64th ed. Available at: http://www.bp.com/content/dam/bp/pdf/Energy-eco-nomics/statistical-review-2015/bp-statistical-review-of-world-energy-2015-full-report.pdf (accessed 03 March 2015).
4. Inaba S. A Two-Phase Code for Protoplanetary Disks. A&A, 2005, vol. 431, pp. 365-379.
5. Hubber D.A. SEREN - a New SPH Code for Star and Planet Formation Simulations. Algorithms and Tests. A&A, 2011, vol. 529. Available at: http://www.aanda.org/articles/aa/pdf/2011/05/ aa14949-10.pdf (accessed 03 March 2015).
6. Vshivkov V.A. Chislennoe reshenie trekhmernykh zadach dina-miki samogravitiruyushchikh mnogofaznykh sistem [Numerical
solution of three-dimensional problems of the dynamics of self-gravitating multiphase systems]. Scientific Bulletin of NSTU, 2011, vol. 44, no. 3, pp. б9-80.
7. Belotserkovskii O.M., Davydov Iu.M. Metod krupnykh chastits v gazovoy dinamike [Method of large particles in gas dynamics]. Moscow, Nauka Publ., 1982. 392 p.
8. Harlow F.H. The Particle-in-Cell Method for Numerical Solution of Problems in Fluid Dynamics. Proc. of Symposium in Applied Mathematics. New York, 1963. Vol. 15, no. 10, pp. 269-288.
9. Kulikovskii A.G., Pogorelov N.V., Semenov A.Yu. Matema-ticheskie voprosy chislennogo resheniya giperbolicheskikh sistem uravneniy [Mathematical problems of numerical solutions of hyperbolic systems of equations]. Moscow, Fizmatlit Publ., 2001.
608 p.
10. Hockney R.W., Eastwood J.W. Computer Simulation Using Particles. New York, Adam Hilger, iOp Publication Ltd., 1989. 540 p.
11. Moscardini L. Cosmological Particle-Mesh N-Body Simulations and the Role of Initial Conditions. Rend. Sem. Mat. Univ. Pol. Torino, 1993, vol. 51, no. 3, pp. 249-266.
69
Протасов В.А., Куликов И.М. PADME - новый код для моделирования процесса формирования георесурсов ... С. 61-70
12. Vshivkov V.A. Computational Methods for Ill-Posed Problems of Gravitational Gasodynamics. Journal of Inverse and Ill-Posed Problems, 2011, vol. 19, no 1, pp. 151-166.
13. Boronina M.A., Vshivkov V.A., Dudnikova G.I. Neyavnaya skhe-ma dlya resheniya uravneniy Maksvella v oblastyakh s razlichny-mi masshtabami [An implicit scheme for solving Maxweell’s equations in domains with different scales]. Reports of the Academy of Sciences of Higher School of the Russian Federation. 2014, no. 4, pp. 39-46.
14. Birdsall C.K., Fuss D. Clouds-in-Clouds, Clouds-in-Cells Physics for Many-Body Plasma Simulation. Journal of Computational Physics, 1997, vol. 135, pp. 141-148.
15. Vshivkov V.A. Parallelnaya realizatsiya na superEVM modeli ga-zovoy komponenty samogravitiruyushchego protoplanetnogo dis-ka [Parallel implementation of the model of gas component of self-gravitating protoplanetary disk on supercomputer]. Computational technologies. 2007, vol. 12, no. 3, pp. 38-52.
16. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Berlin, Heidelberg; New York, Springer-Verlag, 1999. 640 p.
17. Rybakin B.P., Shider N.I. Postroenie parallelnykh algoritmov dlya resheniya zadach gravitatsionnoy gazovoy dinamiki [Development of parallel algorithms for solving the problems of gravitational gas dynamics]. Numerical methods and programming, 2010, vol. 11, pp. 388-394.
18. Rybakin B.P. Chislennoe modelirovanie zadach gazovoy dinamiki na mnogoprotsessornykh EVM [Numerical modelling of gas dynamics problems on multiprocessing computers]. Parallelnye vychi-sleniya i zadachi upravleniya. Doklady pyatoy Mezhdunarodnoy konferentsii [Proc. of the fifth Intern. conf. Parallel computing and control problems]. Moscow, 2010. pp. 161-170.
19. Boreskov A.V., Kharlamov A.A. Osnovy raboty s tekhnologiey CUDA [Basics of CUDA technology]. Moscow, DMK Press, 2010. 232 p.
20. Sanders J., Kandrot E. CUDA by Example. An introduction to general-purpose GPU programming. Moscow, DMK Press, 2011. 232 p.
21. Kulikov I.M., Chernykh I.G., Glinskiy B.M. AstroPhi: program-mny kompleks dlya modelirovaniya dinamiki astrofizicheskikh obektov na gibridnykh superEVM, osnashchennykh uskoritelya-mi Intel Xeon Phi [AstroPhi: a hydrodynamical code for complex modelling of astrophysical objects dynamics by means of hybrid architecture supercomputers on Intel Xeon Phi base]. Bulletin of the South Ural state university. Series «Computational mathematics and software engineering», 2013, vol. 2, no. 4, pp. 57-79.
22. Kulikov I.M. Parallelnaya realizatsiya chislennoy modeli stolkno-veniya galaktik [Parallel program for numerical model of colliding galaxies]. Scientific Bulletin of NSTU, 2011, vol. 9, no. 4, pp. 71-78.
Received: 11 March 2015.
70
