\documentclass{cmi}
\usepackage{upgreek}

\begin{document}


\classify{УДК 519.688}

\renewcommand*{\thefootnote}{\fnsymbol{footnote}}

\title{ОСОБЕННОСТИ РЕАЛИЗАЦИИ АЛГОРИТМА \\ TREECODE ДЛЯ РЕШЕНИЯ ЗАДАЧИ N-ТЕЛ \\ С ИСПОЛЬЗОВАНИЕМ  ГРАФИЧЕСКИХ УСКОРИТЕЛЕЙ\footnote{Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии (ПаВТ) 2021»}%
}

\author{\copyright~2021~г. А.В.~Титов, А.В.~Хоперсков}


\address{Волгоградский Государственный Университет\\
	(400062 Волгоград, пр.~Университетский, д.~100)}

\email{E-mail: \href{mailto:alexandr.titov@volsu.ru}{alexandr.titov@volsu.ru},
\href{mailto:khoperskov@volsu.ru}{khoperskov@volsu.ru}}

\received{Поступила в редакцию: 12.02.2021}

\maketitle{}

\renewcommand*{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}


\begin{abstract}%

Иерархические методы вычисления гравитационных сил для систем $N$-тел позволяют существенно увеличить качество численного моделирования при решении различных астрофизических задач за счет увеличения числа элементов $N$, поскольку вместо вычислительной сложности $\sim O(N^2)$ для прямого метода, мы имеем $N\, \log(N)$ при использовании приближенного метода TreeCode, что позволяет существенно увеличить число частиц в численных моделях. 
Разработано новое программное обеспечение для решения динамической задачи с большим числом частиц для моделирования галактических бесстолкновительных компонент, в частности, звездной подсистемы и темной массы. 
 В работе представлены результаты тестирования алгоритма TreeCode для параллельной реализациии на графических ускорителях NVidia Tesla. Для построения иерархической системы сеток нами реализован быстрый алгоритм построения октодеревьев, основанный на пространственной кривой Мортона.
Для оценок качества построенной численной модели используем для сравнения результаты моделирования на основе прямого вычисления сил взаимодействия между всеми $N$ частицами системы.
Проведен анализ быстродействия различных реализаций алгоритмов решения задачи N-тел и выполнения интегральных законов сохранения физических характеристик для гравитирующих систем. В частности, проанализированы законы сохранения энергии и момента импульса для вращающегося самогравитирующего диска. 
Рассмотрены модели с различными критериями оценки удаленности частицы и значениями угла раскрытия~$\theta$.  


\keywords{задача N-тел,
метод Treecode,
параллельные вычисления,
графические ускорители}
\end{abstract}
\flushleft{\textbf{ОБРАЗЕЦ ЦИТИРОВАНИЯ}}
\justifying
\begin{citationplace}
Титов~А.В., Хоперсков~А.В. Особенности реализации алгоритма Treecode для решения задачи N-тел с использованием графических ускорителей // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2021. Т.~X, №~Y. С.~Z1--Z2.
DOI:~10.14529/cmseXXXXXX.
\end{citationplace}

\markboth{А.В.~Титов, А.В.~Хоперсков}{Особенности реализации алгоритма Treecode для решения задачи N-тел с использованием графических ускорителей \ldots}

\section*{Введение}
\label{sec-intro}

Задача $N$-тел восходит к небесной механике, где решаются  уравнения движения с учетом гравитационного взаимодействия для системы небесных тел \cite{Aarseth-2008book, Greenspan-2004book}. 
С развитием вычислительной техники стало возможным численное интегрирование уравнений движения первоначально для сотен частиц, что позволило моделировать коллективные процессы в приложении к галактическим системам \cite{Ostriker-Peebles-1973bar}, а с переходом к $N\sim 10^6$ объектом исследований стала метагалактика и первоначально в таких космологических моделях под частицей понималась галактика (см., например, \cite{Klypin-etal-1993cosmol}).
Современные астрофизические модели $N$-тел могут содержать свыше 1~млрд элементов (частиц различных видов, моделирующих звезды, темное вещество, газ) \cite{Schaye-2015EAGLE, Pillepich-etal-2018IllustrisTNG}, и даже превышают $10^{10}$ элементов в проекте IllustrisTNG.
Такие проекты, как Illustris, IllustrisTNG, EAGLE (Evolution and Assembly of GaLaxies and their Environments) \cite{Schaye-2015EAGLE}, AGORA \cite{Kim-etal-2014agora}, FIRE (Feedback In Realistic Environments), Indra, Millennium и Millennium-II Simulations \cite{Boylan-Kolchin-etal-2009Millennium-II} и др. являются своего рода галактическими и космологическими <<фабриками>> для симуляций эволюции Вселенной на различных масштабах.
 Результаты таких проектов, например, TNG50, TNG100,  TNG300, AGORA выступают в качестве источников исходных данных, которые позволяют на основе обработки таких симуляций изучать механизмы формирования галактических структур, их морфологию и кинематику \cite{Zhao-etal-2020IllustrisTNG, khoperskov-2021IllustrisTNG, Ubler-etal-2021TNG50}.
 
Модели $N$-тел активно используются для моделирования широкого круга задач в физике плазмы, физической электронике при описании ионных пучков, динамике разреженных газов \cite{Hockney-Eastwood-1987book, Greenspan-2004book}.
К этому классу моделей относяться методы молекулярной динамики для описания эволюции молекулярных систем с использованием как классических, так и квантовохимических моделей \cite{Nemukhin-etal-2009rev, Kholmurodov-etal-2003rev, suplatov-etal-2019mol-gpu}. 

Развитие быстрых методов, в частности, алгоритма Treecode было первоначально связано с моделированием гравитирующих систем \cite{barnes1986}.
Гравитационное взаимодействие $N$-тел является ресурсоемкой задачей из-за дальнодействующего характера взаимодействия и использование прямого метода, имеющего вычислительную сложность $O(N^2)$, не позволяет получить решение за приемлемое время для систем с большим количеством частиц. Алгоритм Treecode позволяет снизить вычислительную сложность до $O(N\log N)$ за счет кластеризации удаленных частиц и использования аппроксимации их влияния на каждую частицу \cite{barnes1986, hamada2009}.

 Важным направлением является развитие технологий программирования гибридных вычислительных систем, которая обеспечивает автоматическую адаптацию кода при изменении вычислительной архитектуры  \cite{Gulenok-etal-2017hybrid}.
Применение комплексного подхода к решению задач моделирования сложных физических систем с эффективным использованием архитектуры суперкомпьютера является насущной задачей \cite{Glinskiy-etal-2017integr}.
Массовый переход на вычислительные системы с графическими процессорами NVIDIA Tesla в последние годы для проведения суперкомпьютерных вычислений ставит задачу реализации эффективных кодов для вычисления парных взаимодействий для системы $N$-тел \cite{Kuksheva-etal-2005paral, Glinskiy-etal-2014paral, suplatov-etal-2019mol-gpu}.

В работе описана реализация параллельного кода для графических ускорителей для численного решения динамической задачи для $N$ гравитационно взаимодействующих частиц, что представляет интерес для решения задач, прежде всего, физики галактик. В основе реализованного Treecode-алгоритма лежит построение дерева разрядов посредством преобразования координат частиц к целочисленным индексам с использованием кривой Мортон, описанное в разделе 1. В разделе~2 проведен анализ производительности программной реализации и рассчитаны изменения со временем углового момента и полной энергии системы $N$-тел. Показана сходимость решений при уменьшении угла раскрытия~$\theta$. В заключении
приводится краткая сводка результатов, полученных в работе.



\section{Модель $N$-тел и параллельная реализация кода}

Рассмотрим систему из $N$ гравитационно взаимодействующих частиц. Уравнения движения для частицы с индексом $i$ ($i = 1, ..., N$) имеют вид:
\begin{equation}\label{eq-1}
    \frac{d\mathbf{r}_i}{dt} = \mathbf{v}_i,
\end{equation}
\begin{equation}\label{eq-2}
   \frac{d\mathbf{v}_i}{dt} = - G \sum^{N}_{j = 1,\, j \ne i}m_j\,\frac{\mathbf{r}_i - \mathbf{r}_j}{|\mathbf{r}_i - \mathbf{r}_j + \varepsilon|^3} \,,
\end{equation}
где вектора $\mathbf{r}$ и $\mathbf{v}$ определяют положение и скорость частицы соответственно, $m$ --- масса частицы, $\varepsilon \ge 0$ --- радиус гравитационного обрезания. Случай $\varepsilon > 0$ позволяет рассматривать бесстолкновительные системы, когда число реальных объектов $N_*$ существенно превышает число частиц $N$ в численной модели. Например, число звезд в типичной крупной галактике $N_*\sim 10^{11}$, а в численной модели мы ограничены, как правило, $N\sim 10^{7}$. 
При моделировании систем электрически заряженных частиц указанное различие может быть сколь-угодно большим.
Для обеспечения бесстолкновительности системы, в которых сила взаимодействия между ближайшими частицами мала по сравнению с вкладом остальных частиц, в численной модели необходимо задавать конечное значение $\varepsilon > 0$ в зависимости от $N$ и пространственной конфигурации изучаемой системы \cite{fridman-2013book, Smirnov-Sotnikova2017}.

Ключевая идея алгоритма Treecode заключается в 
иерархическом разбиении пространства на древовидную структуру. Такая структура позволяет более точно вычислять силы от близкорасположенных частиц, а вклад от более далеких частиц учитывается приближенно так, чтобы иметь контролируемую погрешность суммарной силы.
При вычислении силы для каждой частицы происходит рекурсивный обход ячеек дерева и на каждом уровне оценивается удаленность ячейки от частицы с помощью критерия принятия решения. Существуют различные критерии, позволяющие проверить, достаточно ли удалена частица для оценки силы. 
В работе используются два различных критерия принятия решения: стандартный критерий, введенный в работе~\cite{barnes1986}:
\begin{equation}\label{eq:1}
    d \ge \frac{l}{\theta},
\end{equation}
и модифицированный, который дает более точную оценку силы в случаях, когда центр масс ячейки сильно смещен от ее геометрического центра~\cite{lukat2016}:
\begin{equation}\label{eq:2}
    d \ge \frac{l}{\theta} + \delta,
\end{equation}
где $d$ --- расстояние от частицы до центра масс ячейки, $l$ --- размер ячейки, $\theta$ --- угол раскрытия (opening angle), $\delta$ --- расстояние от центра масс ячейки до ее геометрического центра. Параметр $\theta$ определяет значение максимального угла, под которым мы видим ячейку из точки, для которой вычисляется гравитационная сила.

\fig{height=0.9\textheight}{figTitov_00}{Диаграмма потоков параллельного приложения}


На рис. \ref{figTitov_00} представлена диаграмма потоков разработанного приложения.  Параллельный алгоритм включает в себя следующие этапы:

\begin{itemize}
    \item[(1)] На первом этапе происходит построение дерева разрядов. Листовые узлы такого дерева содержат упорядоченные по возрастанию индексы частиц, а внутренние узлы разделяют, содержащиеся в нем индексы на два линейных диапазона в соотвествии с их общими префиксами в бинарном представлении индексов. Перед построением дерева происходит определение ограничевающего объема посредством вычисления максимального расстояния $R_{\max}$ от центра системы до дальней частицы с использованием алгоритма параллельной редукции. Для всех частиц системы производится преобразование координат частиц к целочисленным значениям, которые используются для вычисления индекса частицы. Индексы частиц определяются на основе пространственной кривой Мортона.
    Используя двоичное представление кодов Мортона строим дерево разрядов, являющееся префиксным деревом, в котором внутренние узлы хранят наибольший общий префикс между ключами в соответствующем поддереве~\cite{karras2010}. 
    
    \item[(2)] На втором этапе происходит построение октодерева. Используя информацию о длинах префиксов в узлах дерева разрядов определяется количество элементов в октодереве~\cite{karras2010}. Далее производится поиск родительских узлов октодерева с использованием ссылки  на предков каждого узла дерева разрядов. При привязке частицы к узлу октодерева эта частица передает ячейке свои характеристики: массу и коррдинаты. Проходя дерево от его листовых узлов к внутренним производим суммирование массы и вычисление центров масс узлов, выполняя синхронизацию потоков после перехода на новый уровень глубины октодерева.
    
    \item[(3)] Ядро для вычисления сил использует нерекурсивный алгоритм обхода дерева. При проходе узлов октодерева с использованием одного из критериев принятия решения (\ref{eq:1}) или (\ref{eq:2}) определятся необходимость использования аппроксимации. 
    
    \item[(4)] Ядро, используемое для интегрирования системы уравнений (\ref{eq-1}), (\ref{eq-2}) включает в себя предиктор и корректор. На этапе предиктор рассчитывается скорость и положение частицы по формулам:
    \begin{equation}\label{eq-5}
        v_{n + 1/2} = v_n + \frac{dv_n}{dt}\frac{dt}{2} \,,
    \end{equation}
    \begin{equation}\label{eq-6}
        x_{n + 1} = x_n + v_{n + 1/2}dt \,,
    \end{equation}
    где $dt$ --- шаг по времени.
Затем выполняется перестроение дерева и вычисляются силы с использованием новых положений частиц. На этапе корректор выполняется уточнение скорости в момент времени $n+1$ по формуле:
     \begin{equation}\label{eq-7}
        v_{n + 1} = v_{n + 1/2} + \frac{dv_{n + 1/2}}{dt}\frac{dt}{2} \,.
    \end{equation}

\end{itemize}

Все тестовые расчеты проводились для бесстолкновительного гравитирующего трехмерного вращающегося диска. Такой диск является равновесным, находясь на балансе гравитационных, центробежной и сил <<давления>> (из-за хаотического движения частиц) в радиальном направлении и на балансе гравитационных и сил <<давления>> в вертикальном направлении (см. подробное описание в \cite{khoperskov-2010, fridman-2013book}).
Хаотическое движение частиц в таких бесстолкновительных системах можно характеризовать тремя пространственными распределениями дисперсии скоростей ($c_r(r), c_\varphi(r), c_z(r)$), соответственно, в радиальном, азимутальном и вертикальном направлениях. 
Мы использовали радиальные профили параметров бесстолкновительного диска из работы \cite{khoperskov2003}, пересчитанные в отсутствии темной массы. 
За счет варьирования начального профиля дисперсии скоростей частиц $c_r(r)$ можно строить как гравитационно устойчивые системы, сохраняющие в целом осесимметричное распределение параметров, так и диски, допускающие развитие гравитационной неустойчивости \cite{fridman-2013book}. В последнем случае в процессе эволюции формируются спиральные структуры значительной амплитуды, что является хорошим тестом для анализа выполнения физических законов сохранения.  

\section{Анализ производительности и законов сохранения}

Параллельные расчеты проводились с использованием графического сопроцессора Nvidia Tesla V100, а тестирование последовательной реализация Treecode~\cite{treecode} проводилось на Intel Xeon E5-2650 v3. Для оценки эффективности распараллеливания проводились сравнения времени вычисления для различных реализаций решения задачи $N$-тел:  параллельная реализация прямого метода (каждая частица с каждой) с использованием программного обеспечения для CUDA~\cite{khrapov2018},  открытая версия алгоритма Treecode от Joshua E. Barnes~\cite{treecode}.  

Были проведены серии расчетов с различным числом частиц и для двух критериев принятия решения. На рис. ~\ref{figTitov_01} показаны графики сравнения времени, необходимой для одной итерации вычисления силы для различного числа частиц и различных реализаций решения задачи $N$-тел: $t_{tc} GPU$ --- время вычислений с использованием разработанного параллельного алгоритма для GPUs, $t_d GPU$ --- время вычислений с использованием прямого алгоритма на GPUs \cite{khrapov2018}, $t_{tc} CPU$ --- время вычислений для последовательной версии Treecode. 

\fig{width=\textwidth}{figTitov_01}{Сравнение времени расчета силы для различных реализаций с теоретическими \\ оценками (слева) и полученное ускорение $S$ (справа)}


На рис. \ref{figTitov_01} справа представлены результаты ускорения вычислений для модели дисковой галактики \cite{khoperskov2003}. Параллельная реализация Treecode показала ускорение вычислений в 40 -- 45 раз по сравнению с последовательной версией, а также ускорение до 68 раз по сравнению с параллельной реализаций прямого алгоритма при использовании NVidia Tesla V100.
Мы провели сравнение точности выполнения основных законов сохранения для выбранной модели диска: углового момента и полной энергии. На рис. \ref{figTitov_02} и \ref{figTitov_03} представлены графики изменения относительного отклонения углового момента (слева) и энергии (справа) для различных значений параметра $\theta$. 

\fig{width=\textwidth}{figTitov_02}{Изменения относительного отклонения углового момента и энергии для \\ вращающегося диска из $N=2^{16}$ частиц и стандартного критерия принятия решения}


Относительная погрешность вычисления углового момента варьируется в диапазоне от 0.0001 до 0.008 в зависимости от $\theta$ при использовании стандартного критерия, в то время как для модифицированного критерия погрешность варьируется в диапазоне от 0.00002 до 0.002. Относительная погрешность вычисления полной энергии для стандартного критерия лежит в диапазоне от 0.0002 до 0.01, а для модифицированного критерия от 0.0001 до 0.003. 
\fig{width=\textwidth}{figTitov_03}{Изменения относительного отклонения углового момента и энергии для диска с $N=2^{16}$ при использовании $\delta$ критерия}


\section*{Заключение}

Методы молекулярной динамики широко используются для моделирования разнообразных физико-технических объектов и опираются на вычисление парных взаимодействий между частицами \cite{Hockney-Eastwood-1987book, Aarseth-2003, Greenspan-2004book, Aarseth-2008book,  Kireyeu-etal-2020-Nuclear-Clusters-Hypernuclei}.
Отметим широкий перечень областей применимости моделей типа $N$-тел, включая физическую кинетику, физику плазмы и ионных пучков, космологию, физику галактик, позволяя моделировать как барионное вещество, так и темную материю в приложении к астрофизическим проблемам \cite{Klypin-etal-1993cosmol, zasov-2008, zasov-2017}.

Наш модуль для вычисления суммарной гравитационной силы от большого числа частиц на основе TreeCode можно использовать для моделирования самых различных галактических компонент, включая динамику гравитирущего газа как на больших масштабах ($\sim 10-10^4$ пк) \cite{Schaye-2015EAGLE, Pillepich-etal-2018IllustrisTNG, Boylan-Kolchin-etal-2009Millennium-II, fridman-2013book}, так и на малых масштабах ($\sim 0.001-10$ пк) \cite{korolev-etal-2020super}. 

В работе представлены предварительные результаты тестирования разработанного параллельного приложения для моделирования динамики $N$ гравитационно взаимодействующих частиц с использованием графических ускорителей Nvidia Tesla. 
В качестве модели для тестирования использовался вращающийся бесстолкновительный самогравитирующийся трехмерный диск, первоначальное состояние которого допускало развитие гравитационной неустойчивости, на нелинейной стадии которой формировались центральная перемычка (так называемая бар-мода) и спиральный узор. 
Разработанная параллельная реализация алгоритма Treecode для GPUs обеспечивает приемлемое выполнение интегральных законов сохранения энергии и углового момента в самогравитирующем диске на уровне, сопоставимом с расчетами для CPU. 

Производительность программного обеспечения замерялась на графическом сопроцессоре NVidia Tesla V100 (Volta). 
В условиях проведения тестов величина ускорение составила 45 раз по сравнению с реализацией~\cite{treecode} и до 68 раз по сравнению с параллельным алгоритмом на основе прямого вычисления гравитационной силы~\cite{khrapov2018}. 


\vspace{1em}
{\it Работа выполнена в рамках государственного задания Министерства науки и высшего образования Российской Федерации №0633-2020-0003. Расчеты проводились на оборудовании ЦКП <<Суперкомпьютерный центр коллективного пользования ВолГУ>>.}

\begin{biblio}

\bibitem{Aarseth-2008book}
Aarseth S.J. The Cambridge N-body lectures. Lecture Notes in Physics. 2008. Vol. 760. 402~p.
DOI:~\href{http://dx.doi.org/10.1007/978-1-4020-8431-7}{10.1007/978-1-4020-8431-7}.

\bibitem{Greenspan-2004book}
Greenspan D. N-body problems and models. World Scientific, 2004. 193~p.

\bibitem{Ostriker-Peebles-1973bar}
Ostriker~J.P., Peebles~P.J.E. A Numerical Study of the Stability of Flattened Galaxies: or, can Cold Galaxies Survive? // Astrophysical Journal. 1973. Vol.~186. P.~467--480. 
DOI:~\href{http://dx.doi.org/10.1086/152513}{10.1086/152513}.

\bibitem{Klypin-etal-1993cosmol}
Klypin A., Holtzman J., Primack J., Regos E. Structure Formation with Cold plus Hot Dark Matter // Astrophysical Journal. 1993. Vol.~416. P.~1--16. 
DOI:~\href{https://doi.org/10.1086/173210}{10.1086/173210}.

\bibitem{Schaye-2015EAGLE}
Schaye J., Crain R.A., Bower R.G., et al.  The EAGLE project: simulating the evolution and assembly of galaxies and their environments // Monthly Notices of the Royal Astronomical Society. 2014. Vol.~446. P.~521--554. 
DOI:~\href{https://doi.org/10.1093/mnras/stu2058}{10.1093/mnras/stu2058}.


\bibitem{Pillepich-etal-2018IllustrisTNG}
Pillepich A., Springel V., Nelson D., et al. Simulating galaxy formation with the IllustrisTNG model // Monthly Notices of the Royal Astronomical Society. 2018. Vol.~473. P.~4077--4106.
DOI:~\href{https://doi.org/10.1093/mnras/stx2656}{10.1093/mnras/stx2656}.

\bibitem{Kim-etal-2014agora}
Kim J., Abel T., Agertz O., et al. The AGORA High-resolution Galaxy Simulations Comparison Project // The Astrophysical Journal Supplement. 2014. Vol.~210.  id. 14. 20 P. 
DOI:~\href{https://doi.org/10.1088/0067-0049/210/1/14}{10.1088/0067-0049/210/1/14}.

\bibitem{Boylan-Kolchin-etal-2009Millennium-II}
Boylan-Kolchin M., Springel V., White S.D.M., Jenkins A., Lemson G. Resolving cosmic structure formation with the Millennium-II Simulation // Monthly Notices of the Royal Astronomical Society. 2009. Vol.~398. P.~1150--1164. 
DOI:~\href{https://doi.org/110.1111/j.1365-2966.2009.15191.x}{10.1111/j.1365-2966.2009.15191.x}.

\bibitem{Zhao-etal-2020IllustrisTNG}
Zhao D., Du M., Ho, L.C., Debattista V.P., Shi J.
Barred Galaxies in the IllustrisTNG Simulation // The Astrophysical Journal. 2020. Vol.~904. id.~170. 12~P.
DOI:~\href{https://doi.org/10.3847/1538-4357/abbe1b }{10.1111/10.3847/1538-4357/abbe1b}.

\bibitem{khoperskov-2021IllustrisTNG}
Khoperskov S., Zinchenko I., Avramov B. et al.
Extreme kinematic misalignment in IllustrisTNG galaxies: the origin, structure, and internal dynamics of galaxies with a large-scale counterrotation // Monthly Notices of the Royal Astronomical Society. 2021. Vol.~500.  P.~3870--3888. DOI:~\href{http://dx.doi.org/10.1093/mnras/staa3330}{10.1093/mnras/staa3330}.  

\bibitem{Ubler-etal-2021TNG50}
Ubler H., Genel S., Sternberg A., Genzel R. et al. The kinematics and dark matter fractions of TNG50 galaxies at $z = 2$ from an observational perspective // Monthly Notices of the Royal Astronomical Society. 2021. Vol.~500.  P.~4597--4619.  
DOI:~\href{http://dx.doi.org/10.1093/mnras/staa3464}{10.1093/mnras/staa3464}.  

\bibitem{Hockney-Eastwood-1987book}
%Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987. 640 с. 
Hockney R.W., Eastwood J.W. Computer Simulation Using Particles. CRC Press, 1988, 540~p. ISBN 9780852743928.

\bibitem{Nemukhin-etal-2009rev}
Nemukhin A.V., Grigorenko B.L., Savitsky A.P.
Computer Modeling of the Structure and Spectra of Fluorescent Proteins // Acta Naturae. 2009. Vol.~1, no.~2. P.~33--43.
DOI:~\href{https://doi.org/10.32607/20758251-2009-1-2-33-43}{10.32607/20758251-2009-1-2-33-43}.  

\bibitem{Kholmurodov-etal-2003rev}
Холмуродов Х.Т., Алтайский М.В., Пузынин И.В., Дардин Т., Филатов Ф.П. Методы молекулярной динамики для моделирования физических и биологических процессов // Физика элементарных частиц и атомного ядра. 2003. Т. 34. № 2. С. 474--515
  
\bibitem{suplatov-etal-2019mol-gpu}
Суплатов Д.А., Шарапова Я.А., Попова Н.Н., Копылов К.Е., Воеводин Вл.В., Швядас В.К. Молекулярная динамика в силовом поле FF14SB в воде TIP4P-Ew, и в силовом поле FF15IPQ в воде SPC/Eb: сравнительный анализ на GPU и CPU // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2019. Т. 8, № 1. С. 71–88.
DOI:~\href{http://dx.doi.org/10.14529/cmse190105}{10.14529/cmse190105}.  

\bibitem{barnes1986} Barnes~J., Hut~P. A hierarchical O(N log N) force-calculation algorithm // Nature. 1986. Vol.~324. P.~446--449. 
DOI:~\href{https://doi.org/10.1038/324446a0}{10.1038/324446a0}.

\bibitem{hamada2009} Hamada~T. Nitadori~K., Benkrid K., et al. A novel multiple-walk parallel algorithm for the Barnes–Hut treecode on GPUs -- towards cost effective, high performance N-body simulation // Computer Science -- Research and Development volume. 2009. Vol.~24. P.~21--31. 
DOI:~\href{https://doi.org/10.1007/s00450-009-0089-1}{10.1007/s00450-009-0089-1}.

\bibitem{Gulenok-etal-2017hybrid}
Gulenok A.A., Dordopulo A.I., Levin I.I., Gudkov V.A. Hybrid computer system programming technology with adaptation and scaling of calculations. Bulletin of the South Ural State
University. Series: Computational Mathematics and Software Engineering. 2017. Vol.~6, no.~1. P.~73--86. 
DOI:~\href{https://doi.org/10.14529/cmse170105}{10.14529/cmse170105}.

\bibitem{Glinskiy-etal-2017integr}
Glinskiy B., Kulikov I., Chernykh I., Snytnikov A., Sapetina A., Weins D. The integrated approach to solving large-size physical problems on supercomputers // Communications in Computer and Information Science. 2017. Vol.~793. P.~278--289. 
DOI:~\href{https://doi.org/10.1007/978-3-319-71255-0\_22}{10.1007/978-3-319-71255-0\_22}.

\bibitem{Kuksheva-etal-2005paral}
Kuksheva, E.A., Malyshkin V.E., Nikitin S.A., Snytnikov A.V., Snytnikov V.N., Vshivkov V.A. Supercomputer simulation of self-gravitating media // Future Generation Computer Systems. 2005. Vol.~21. P.~749--757.
DOI:~\href{https://doi.org/10.1016/j.future.2004.05.019}{10.1016/j.future.2004.05.019}.

\bibitem{Glinskiy-etal-2014paral}
Glinskiy B.M., Kulikov I.M., Snytnikov A.V., Romanenko A.A., Chernykh I.G., Vshivkov V.A. Co-design of parallel numerical methods for plasma physics and astrophysics // Supercomputing Frontiers and Innovations. 2014. Vol.~1. P.~88--98. 
DOI:~\href{https://doi.org/10.14529/jsfi140305}{10.14529/jsfi140305}.

\bibitem{fridman-2013book}
Fridman A.M., Khoperskov A.V. Physics of Galactic Disks. Cambridge International Science Publishing Ltd, 2013. 754 p.

\bibitem{Smirnov-Sotnikova2017}
Smirnov, A.A., Sotnikova, N.Ya., Koshkin, A.A. Simulations of slow bars in anisotropic disk systems // Astronomy Letters. 2017. Vol.~43 P.~61--74. 
DOI:~\href{https://doi.org/10.1134/S1063773717020062}{10.1134/S1063773717020062}.

\bibitem{lukat2016} Lukat G., Banerjee R. A GPU accelerated Barnes–Hut tree code for FLASH4 // New Astronomy. 2016. Vol.~45. P.~14--28. 
DOI:~\href{https://doi.org/10.1016/j.newast.2015.10.007}{10.1016/j.newast.2015.10.007}.

\bibitem{karras2010} Karras~T. Maximizing Parallelism in the Construction of BVHs, Octrees, and K-d Trees // Proceedings of the 4th Symposium on High-Performance Graphics.  2012.  P.~33--37. DOI:~\href{http://dx.doi.org/10.2312/EGGH/HPG12/033-037}{10.2312/EGGH/HPG12/033-037}.

\bibitem{khoperskov-2010}
Khoperskov A.V., Bizyaev D., Tiurina N., Butenko M. Numerical modelling of the vertical structure and dark halo parameters in disc galaxies // Astronomische Nachrichten. 2010. Vol.~331. P.~731-745.
DOI:~\href{http://dx.doi.org/10.1002/asna.200911402}{10.1002/asna.200911402}.

\bibitem{khoperskov2003} Khoperskov~A.V., Tyurina~N.V. A Dynamical Model of the Galaxy // Astronomy Reports. 2003.  Vol.~47. P.~443--457.
DOI:~\href{https://doi.org/10.1134/1.1583771}{10.1134/1.1583771}.

\bibitem{khrapov2018} Khrapov~S.S., Khoperskov~S.A., Khoperskov~A.V. New features of parallel implementation of N-body problems on GPU // Bulletin of the South Ural State University, Series: Mathematical Modelling, Programming and Computer Software. 2018. Vol.~11. No.~1. P.~124--136.
DOI:~\href{https://dx.doi.org/10.14529/mmp180111}{10.14529/mmp180111}.

\bibitem{treecode}
Treecode guide. URL:~\href{https://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html}{https://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html} (дата обращения: 30.11.2020).

\bibitem{Aarseth-2003}
Aarseth S.J. Gravitational N-Body Simulations Tools and Algorithms. Cambridge University Press, 2003. 431~p.

\bibitem{Kireyeu-etal-2020-Nuclear-Clusters-Hypernuclei}
Kireyeu V., Aichelin J., Bratkovskaya, E., Le Fevre A., Lenivenko V., Kolesnikov V., Leifels Y., Voronyuk V. PHQMD Model for the Formation of Nuclear Clusters and Hypernuclei in Heavy Ion // 
Bulletin of the Russian Academy of Sciences: Physics. 2020. Vol.~84. P.~957--961. DOI:~\href{https://doi.org/10.3103/S1062873820080171}{10.3103/S1062873820080171}.

\bibitem{zasov-2008}
Zasov A.V., Moiseev A.V., Khoperskov A.V., Sidorova E.A. Early-type disk galaxies: structure and kinematics // Astronomy Reports. 2008. Vol.~52. P.~79--93. 
DOI:~\href{https://doi.org/10.1134/S1063772908020017}{10.1134/S1063772908020017}.
 
\bibitem{zasov-2017}
Zasov A.V., Saburova A.S., Khoperskov A.V., Khoperskov S.A. Dark matter in galaxies // Physics-Uspekhi. 2017. Vol.~60, no.1.  P.~3--39. 
DOI:~\href{https://doi.org/10.3367/UFNe.2016.03.037751}{10.3367/UFNe.2016.03.037751}.

\bibitem{korolev-etal-2020super}
Korolev V.V., Eremin M.A., Kovalenko I.G., Zankovich A.M. Numerical simulations of instability in the shell of a supernova remnant expanding in a weakly inhomogeneous interstellar medium // Mathematical Physics and Computer Simulation. 2020. Vol.~23 No.~3. P.~23---35.
DOI:~\href{https://doi.org/10.15688/mpcm.jvolsu.2020.3.3}{10.15688/mpcm.jvolsu.2020.3.3}.
\end{biblio}


Александр Викторович Титов,  аспирант, кафедра информационных систем и компьютерного моделирования, Волгоградский государственный университет (Волгоград, Российская Федерация), e-mail:~alexandr.titov@volsu.ru

Александр Валентинович Хоперсков, д.ф.-м.н., профессор, кафедра информационных систем и компьютерного моделирования, Волгоградский государственный университет (Волгоград, Российская Федерация), e-mail:~khoperskov@volsu.ru



{\phantom{eee} \hrule \vskip 4 pt}
\classify{} % Код MSC не нужен
%\classify{MSC 00000} % Коды MSC см. http://www.ams.org/msc/

\title{IMPLEMENTATION FEATURES OF THE TREECODE ALGORITHM FOR SOLVING N-BODY PROBLEMS ON GPUs}

\author{\copyright~2019 A.V.~Titov, A.V.~Khoperskov}

\address{Volgograd State University (pr. Universitetsky 100, Vologograd, 400062 Russia)}  

\email{E-mail: \href{mailto:alexandr.titov@volsu.ru}{alexandr.titov@volsu.ru},
\href{mailto:khoperskov@volsu.ru}{khoperskov@volsu.ru}}

\received{Received: 12.02.2021}

\maketitle{}

\begin{abstract}%
Hierarchical methods for calculating gravitational forces in a $N$-body system significantly increase the quality of numerical simulations when solving various astrophysical problems by increasing the number of $N$ elements, since we have the computational complexity $N \, \log(N)$ for the TreeCode approximate method instead of $\sim O(N^2)$ for the direct method, which allows to greatly increase the number of particles in the models. 
We developed new software for solving a dynamic problem with a large number of particles for modeling the collisionless components of the galaxies, in particular, stellar subsystem and dark matter.
The paper presents the test results for the parallel implementation of the Treecode algorithm for the NVidia Tesla GPUs.
To construct a hierarchical grid structure, we implemented a fast parallel octree-construction algorithm based on Morton's space-filling curve.
To assess the quality of the constructed numerical model we use the simulation results based on the direct calculation of the interaction forces between all $N$ particles of the system.
We have compared the performance of the different implementations of algorithms for solving the N-body problem and an analisis  of the fulfillment of the integral physical conservation laws of a self-gravitational system. 
The analysis of the fulfillment of the conservation laws of total energy and angular momentum is carried out for a rotating self-gravitating disk. 
Models with different criteria for a particle remoteness and value of the opening angle $\theta$ are considered.

	\keywordsen{N-body problem, treecode, parallel computing, GPUs}
\end{abstract}

\begin{flushleft}
	\textbf{FOR CITATION}
\end{flushleft}
\justifying

\begin{citationplace}
Titov~A.V., Khoperskov~A.V. Features of the implementation of the TreeCode algorithm for solving N-body problems on GPUs. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2021. Vol.~X, no.~Y. P.~Z1--Z2. (in Russian) DOI: 10.14529/cmseXXXXXX.
\end{citationplace}



\begin{biblio_lat}
%\itemsep -2pt

\bibitem{Aarseth-2008book}
Aarseth S.J. The Cambridge N-body lectures. Lecture Notes in Physics. 2008. Vol. 760. 402 P.
DOI:~\href{http://dx.doi.org/10.1007/978-1-4020-8431-7}{10.1007/978-1-4020-8431-7}.

\bibitem{Greenspan-2004book}
Greenspan D. N-body problems and models. World Scientific, 2004. 193~p.

\bibitem{Ostriker-Peebles-1973bar}
Ostriker~J.P., Peebles~P.J.E. A Numerical Study of the Stability of Flattened Galaxies: or, can Cold Galaxies Survive? // Astrophysical Journal. 1973. Vol.~186. P.~467--480. 
DOI:~\href{http://dx.doi.org/10.1086/152513}{10.1086/152513}.

\bibitem{Klypin-etal-1993cosmol}
Klypin A., Holtzman J., Primack J., Regos E. Structure Formation with Cold plus Hot Dark Matter // Astrophysical Journal. 1993. Vol.~416. P.~1--16. 
DOI:~\href{https://doi.org/10.1086/173210}{10.1086/173210}.

\bibitem{Schaye-2015EAGLE}
Schaye J., Crain R.A., Bower R.G., et al.  The EAGLE project: simulating the evolution and assembly of galaxies and their environments // Monthly Notices of the Royal Astronomical Society. 2014. Vol.~446. P.~521--554. 
DOI:~\href{https://doi.org/10.1093/mnras/stu2058}{10.1093/mnras/stu2058}.


\bibitem{Pillepich-etal-2018IllustrisTNG}
Pillepich A., Springel V., Nelson D., et al. Simulating galaxy formation with the IllustrisTNG model // Monthly Notices of the Royal Astronomical Society. 2018. Vol.~473. P.~4077--4106.
DOI:~\href{https://doi.org/10.1093/mnras/stx2656}{10.1093/mnras/stx2656}.

\bibitem{Kim-etal-2014agora}
Kim J., Abel T., Agertz O., et al. The AGORA High-resolution Galaxy Simulations Comparison Project // The Astrophysical Journal Supplement. 2014. Vol.~210.  id. 14. 20 P. 
DOI:~\href{https://doi.org/10.1088/0067-0049/210/1/14}{10.1088/0067-0049/210/1/14}.

\bibitem{Boylan-Kolchin-etal-2009Millennium-II}
Boylan-Kolchin M., Springel V., White S.D.M., Jenkins A., Lemson G. Resolving cosmic structure formation with the Millennium-II Simulation // Monthly Notices of the Royal Astronomical Society. 2009. Vol.~398. P.~1150--1164. 
DOI:~\href{https://doi.org/110.1111/j.1365-2966.2009.15191.x}{10.1111/j.1365-2966.2009.15191.x}.

\bibitem{Zhao-etal-2020IllustrisTNG}
Zhao D., Du M., Ho, L.C., Debattista V.P., Shi J.
Barred Galaxies in the IllustrisTNG Simulation // The Astrophysical Journal. 2020. Vol.~904. id.~170. 12~P.
DOI:~\href{https://doi.org/10.3847/1538-4357/abbe1b }{10.1111/10.3847/1538-4357/abbe1b}.

\bibitem{khoperskov-2021IllustrisTNG}
Khoperskov S., Zinchenko I., Avramov B. et al.
Extreme kinematic misalignment in IllustrisTNG galaxies: the origin, structure, and internal dynamics of galaxies with a large-scale counterrotation // Monthly Notices of the Royal Astronomical Society. 2021. Vol.~500.  P.~3870--3888. DOI:~\href{http://dx.doi.org/10.1093/mnras/staa3330}{10.1093/mnras/staa3330}.  

\bibitem{Ubler-etal-2021TNG50}
Ubler H., Genel S., Sternberg A., Genzel R. et al. The kinematics and dark matter fractions of TNG50 galaxies at $z = 2$ from an observational perspective // Monthly Notices of the Royal Astronomical Society. 2021. Vol.~500.  P.~4597--4619.  
DOI:~\href{http://dx.doi.org/10.1093/mnras/staa3464}{10.1093/mnras/staa3464}.  

\bibitem{Hockney-Eastwood-1987book}
R.W. Hockney, J.W. Eastwood. Computer simulation using particles. Adam Hilger, Bristol, NY: IOP Publishing Ltd, 1988. 540 P.

\bibitem{Nemukhin-etal-2009rev}
Nemukhin A.V., Grigorenko B.L., Savitsky A.P.
Computer Modeling of the Structure and Spectra of Fluorescent Proteins // Acta Naturae. 2009. Vol.~1, no.~2. P.~33--43.
DOI:~\href{https://doi.org/10.32607/20758251-2009-1-2-33-43}{10.32607/20758251-2009-1-2-33-43}.  

\bibitem{Kholmurodov-etal-2003rev}
Kholmurodov K. T., Altaisky M. V., Puzynin I. V., Darden T., Filatov F. P.
Molecular Dynamics Methods for Simulation of Physical and Biological Processes // Physics of Elementary Particles and Atomic Nuclei. 2003. Vol.~34, no. 2. P. 474--515.  (in Russian) 

\bibitem{suplatov-etal-2019mol-gpu}
Suplatov D.A., Sharapova Y.A., Popova N.N., Kopylov K.E., Voevodin Vl.V., Svedas V.K. Molecular Dynamics in the Force Field FF14SB in Water TIP4P-Ew, and in the Force Field FF15IPQ in Water SPC/Eb: a Comparative Analysis on GPU and CPU // Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2019. Vol.~8, no. 1. P. 71–88. (in Russian) 
DOI:~\href{http://dx.doi.org/10.14529/cmse190105}{10.14529/cmse190105}.  

\bibitem{barnes1986} Barnes~J., Hut~P. A hierarchical O(N log N) force-calculation algorithm // Nature. 1986. Vol.~324. P.~446--449. 
DOI:~\href{https://doi.org/10.1038/324446a0}{10.1038/324446a0}.

\bibitem{hamada2009} Hamada~T. Nitadori~K., Benkrid K., et al. A novel multiple-walk parallel algorithm for the Barnes–Hut treecode on GPUs – towards cost effective, high performance N-body simulation // Computer Science - Research and Development volume. 2009. Vol.~24. P.~21--31. 
DOI:~\href{https://doi.org/10.1007/s00450-009-0089-1}{10.1007/s00450-009-0089-1}.

\bibitem{Gulenok-etal-2017hybrid}
Gulenok A.A., Dordopulo A.I., Levin I.I., Gudkov V.A. Hybrid computer system programming technology with adaptation and scaling of calculations. Bulletin of the South Ural State
University. Series: Computational Mathematics and Software Engineering. 2017. Vol.~6, no.~1. P.~73--86. 
DOI:~\href{https://doi.org/10.14529/cmse170105}{10.14529/cmse170105}.

\bibitem{Glinskiy-etal-2017integr}
Glinskiy B., Kulikov I., Chernykh I., Snytnikov A., Sapetina A., Weins D. The integrated approach to solving large-size physical problems on supercomputers // Communications in Computer and Information Science. 2017. Vol.~793. P.~278--289. 
DOI:~\href{https://doi.org/10.1007/978-3-319-71255-0\_22}{10.1007/978-3-319-71255-0\_22}.

\bibitem{Kuksheva-etal-2005paral}
Kuksheva, E.A., Malyshkin V.E., Nikitin S.A., Snytnikov A.V., Snytnikov V.N., Vshivkov V.A. Supercomputer simulation of self-gravitating media // Future Generation Computer Systems. 2005. Vol.~21. P.~749--757.
DOI:~\href{https://doi.org/10.1016/j.future.2004.05.019}{10.1016/j.future.2004.05.019}.

\bibitem{Glinskiy-etal-2014paral}
Glinskiy B.M., Kulikov I.M., Snytnikov A.V., Romanenko A.A., Chernykh I.G., Vshivkov V.A. Co-design of parallel numerical methods for plasma physics and astrophysics // Supercomputing Frontiers and Innovations. 2014. Vol.~1. P.~88--98. 
DOI:~\href{https://doi.org/10.14529/jsfi140305}{10.14529/jsfi140305}.

\bibitem{fridman-2013book}
Fridman A.M., Khoperskov A.V. Physics of Galactic Disks. Cambridge International Science Publishing Ltd, 2013. 754 p.

\bibitem{Smirnov-Sotnikova2017}
Smirnov, A.A., Sotnikova, N.Ya., Koshkin, A.A. Simulations of slow bars in anisotropic disk systems // Astronomy Letters. 2017. Vol.~43 P.~61--74. 
DOI:~\href{https://doi.org/10.1134/S1063773717020062}{10.1134/S1063773717020062}.

\bibitem{lukat2016} Lukat G., Banerjee R. A GPU accelerated Barnes–Hut tree code for FLASH4 // New Astronomy. 2016. Vol.~45. P.~14--28. 
DOI:~\href{https://doi.org/10.1016/j.newast.2015.10.007}{10.1016/j.newast.2015.10.007}.

\bibitem{karras2010} Karras~T. Maximizing Parallelism in the Construction of BVHs, Octrees, and K-d Trees // Proceedings of the 4th Symposium on High-Performance Graphics.  2012.  P.~33--37. DOI:~\href{http://dx.doi.org/10.2312/EGGH/HPG12/033-037}{10.2312/EGGH/HPG12/033-037}.

\bibitem{khoperskov-2010}
Khoperskov A.V., Bizyaev D., Tiurina N., Butenko M. Numerical modelling of the vertical structure and dark halo parameters in disc galaxies // Astronomische Nachrichten. 2010. Vol.~331. P.~731-745.
DOI:~\href{http://dx.doi.org/10.1002/asna.200911402}{10.1002/asna.200911402}.

\bibitem{khoperskov2003} Khoperskov~A.V., Tyurina~N.V. A Dynamical Model of the Galaxy // Astronomy Reports. 2003.  Vol.~47. P.~443--457.
DOI:~\href{https://doi.org/10.1134/1.1583771}{10.1134/1.1583771}.

\bibitem{khrapov2018} Khrapov~S.S., Khoperskov~S.A., Khoperskov~A.V. New features of parallel implementation of N-body problems on GPU // Bulletin of the South Ural State University, Series: Mathematical Modelling, Programming and Computer Software. 2018. Vol.~11. No.~1. P.~124--136.
DOI:~\href{https://dx.doi.org/10.14529/mmp180111}{10.14529/mmp180111}.

\bibitem{treecode}
Treecode guide. URL:~\href{https://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html}{https://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html} (accessed: 30.11.2020).

\bibitem{Aarseth-2003}
Aarseth S.J. Gravitational N-Body Simulations Tools and Algorithms. Cambridge University Press, 2003. 431~p.

\bibitem{Kireyeu-etal-2020-Nuclear-Clusters-Hypernuclei}
Kireyeu V., Aichelin J., Bratkovskaya, E., Le Fevre A., Lenivenko V., Kolesnikov V., Leifels Y., Voronyuk V. PHQMD Model for the Formation of Nuclear Clusters and Hypernuclei in Heavy Ion // 
Bulletin of the Russian Academy of Sciences: Physics. 2020. Vol.~84. P.~957--961. DOI:~\href{https://doi.org/10.3103/S1062873820080171}{10.3103/S1062873820080171}.

\bibitem{zasov-2008}
Zasov A.V., Moiseev A.V., Khoperskov A.V., Sidorova E.A. Early-type disk galaxies: structure and kinematics // Astronomy Reports. 2008. Vol.~52. P.~79--93. 
DOI:~\href{https://doi.org/10.1134/S1063772908020017}{10.1134/S1063772908020017}.
 
\bibitem{zasov-2017}
Zasov A.V., Saburova A.S., Khoperskov A.V., Khoperskov S.A. Dark matter in galaxies // Physics-Uspekhi (Advances in Physical Sciences). 2017. Vol.~60, no.1.  P.~3--39. 
DOI:~\href{https://doi.org/10.3367/UFNe.2016.03.037751}{10.3367/UFNe.2016.03.037751}.

\bibitem{korolev-etal-2020super}
Korolev V.V., Eremin M.A., Kovalenko I.G., Zankovich A.M. Numerical simulations of instability in the shell of a supernova remnant expanding in a weakly inhomogeneous interstellar medium // Mathematical Physics and Computer Simulation. 2020. Vol.~23 No.~3. P.~23---35.
DOI:~\href{https://doi.org/10.15688/mpcm.jvolsu.2020.3.3}{10.15688/mpcm.jvolsu.2020.3.3}.
\end{biblio_lat}


\end{document}

