CC BY
33
ISSN 1992-6502 (Print)
2016. Т. 20, № 3 (73). С. 137-142
-"Boogie QjrAQnQj
ISSN 2225-2789 (Online)
http://journal.ugatu.ac.ru
УДК 519.63
Масштабируемый параллельный алгоритм для моделирования трехмерной динамики гравитирующих систем
методом частиц
Н. В. Снытников
Институт вычислительной математики и математической геофизики СО РАН Новосибирский государственный университет
Поступила в редакцию 17.03.2016
Аннотация. Предложен параллельный алгоритм для численного решения системы уравнений Власова-Пуассона методом частиц. Используется новый метод динамической балансировки процессоров, которые распределяются между подобластями в соответствии с реальным числом модельных частиц. Метод декомпозиции области позволяет комбинировать сеточный (эйлеров) метод для решения уравнения Пуассона с лагранжевым методом частиц для решения уравнения Власова. Он учитывает физические особенности задач моделирования нестационарных вращающихся дисков (двумерных и трехмерных).
Ключевые слова: гравитирующие системы, уравнение Пуассона, уравнение Власова, метод частиц, динамическая балансировка
Суперкомпьютерное моделирование динамики гравитирующих систем (таких как галактики или околозвездные диски) является одной из наиболее трудных и актуальных проблем в современной вычислительной астрофизике.
Хотя однопроцессорные численные методы были созданы и успешно используются уже более 30 лет [1-5], разработка соответствующих им параллельных алгоритмов требует значительных модификаций, связанных с адаптацией к новым аппаратным и программным архитектурам суперкомпьютеров.
Другие трудности возникают из-за разнообразия моделируемых астрофизических процессов [6-8], для которых существуют специализированные программные пакеты [9-11], не всегда применимые для решения смежных задач.
В данной статье предложен параллельный алгоритм для моделирования вращающихся бес-столкновительных гравитирующих систем методом частиц в ячейках и приведены результаты численных экспериментов (детальное описание алгоритмов приведено в [12]).
Алгоритм может быть применен как для модели тонкого диска (когда отсутствует движение вещества в вертикальном направлении [13]), так и для полностью трехмерной модели [14]. Новым в этом подходе является метод динамической балансировки загрузки процессоров при вычислении траектории модельных частиц, которые в случае вращающихся систем могут многократно переходить между подобластями (и соответствующим им процессорам).
Решаемая система уравнений состоит из уравнения Власова для функции распределения вещества / = /(t, г, и), зависящей от времени t, пространственной координаты г = (х, у, z) и вектора скорости u = (u, v, w), с заданным начальным значением /0(г, и), (в виде некоторого вращающегося диска)
а/ а/ а/
at ar аи
/(0,г,и) = /0(г,и)
Работа поддержана грантами РФФИ №16-07-00916, №16-31-00301, и грантом Президента РФ МК-5915.2016.1.
Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2016».
и уравнения Пуассона для гравитационного потенциала Ф = Ф(^ г) с краевыми условиями для изолированных систем
где G
ДФ(£,г) = 4n:Gp(t, г), Ф(С,г)||гИм = 0, - гравитационная постоянная.
Замыкает систему следующее уравнение для плотности р = р(£, г), которое связывает ее с функцией распределения /:
p(t, г) = J /(t, г, u) du.
Для решения уравнения Власова используется метод частиц в ячейках [1, 2], а для решения уравнения Пуассона - метод свертки. Соответствующий параллельный алгоритм должен подразделять исходную вычислительную область на такие подобласти, что каждая из них могла бы быть обработана одним вычислительным устройством с 1-4 ГБ памяти, и распределять модельные частицы между процессорами так, чтобы они могли иметь доступ к значениям требуемых сеточных функций в подобластях, свободно двигаться между подобластями, и общее количество частиц, передаваемых между процессорами, было минимальным.
В типичном случае максимальное количество частиц, которое можно поместить в одно вычислительное устройство (CPU, GPU, Intel Phi), составляет 8-16 миллионов, а максимальное количество узлов сетки составляет 256 X 256 X 256 для трехмерных задач и 4096 X 4096 для двумерных.
Это означает, что большие вычислительные области (такие как 1024 X 1024 X 1024) должны подразделяться на меньшие подобласти. В то же самое время модельные частицы могут многократно перелетать из одной области в другую в течение одного численного эксперимента (например, если рассматривается моделирование дисков на временных масштабах порядка десятков оборотов, где каждая частица двигается по эллиптической или эпициклической орбите вокруг общего центра масс), что создает проблему пересылок данных на каждом временном шаге. Кроме того, частицы могут быть распределены неравномерно между подобластями.
Например, во вращающемся диске могут возникать гравитационные неустойчивости, которые приводят к появлению разнообразных структур - сгустков вещества или спиральных волн,
которые способны передвигаться по всей вычислительной подобласти. Плотность вещества (и количество частиц) в этих структурах может быть на порядки больше фоновых значений. Следовательно, количество частиц в каждой подобласти может так же отличаться на порядки и меняться со временем, особенно в тех случаях, когда сгусток, содержащий большое количество частиц, двигается вокруг центра масс.
Разработанный параллельный алгоритм использует следующие свойства задачи и численного метода:
• Трудоемкость решения уравнения Власова (т.е. интегрирование траекторий Р1С-ча-стиц) существенно выше (от 10 до 100 раз), чем трудоемкость вычисления гравитационного потенциала. Это связано с тем, что среднее число модельных частиц на одну ячейку сетки в типичных расчетах находится в диапазоне 10-100.
• Чтобы удовлетворять условию Куранта устойчивости численного метода, частица может перелететь за один шаг не далее, чем в соседнюю ячейку. Это означает, что число частиц, которые должны перемещаться между смежными подобластями, намного меньше, чем общее число частиц в смежных подобластях (порядка 10 для 2D задач и 10пуп2 для 3D задач, где пу, п2 -число всех граничных узлов между подобластями).
Общая схема параллельного алгоритма описана в разделе 2, а детальное описание алгоритмов приведено в статье [12]. Раздел 3 посвящен описанию метода вычисления гравитационного потенциала. В разделе 4 приводятся результаты экспериментов.
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ
Исходная вычислительная область (2D или 3D) равномерно подразделяется на К подобластей одинакового размера 51,52, ..., в направлении X. Группа процессоров назначается каждой подобласти (рис. 1). Количество процессоров в группе равно Р^-, Р^- >1. Общее количество процессоров равно Р = £ Р^-. Количество узлов (ячеек) в каждой вычислительной подобласти выбирается таким образом, что все сеточные функции (гравитационный потенциал, плотность, сила) должны полностью помещаться в оперативную память одного вычислительного устройства и оставить в памяти достаточно места для хранения массивов с частицами (т.е. их пространственных координат и скоростей). Обычно это число равно 0,5-2 миллионов узлов (соответствует сетке 1283 или 10242).
U
Количество процессоров в группе выбирается так, чтобы обеспечить приблизительно одинаковое количество частиц на одном процессоре (см. Алгоритм 2 в [12]). Это означает, что если подобласть содержит больше частиц, чем подобласть то число процессоров Р^ в группе будет больше или равно количеству в группе .
В каждой группе процессоров и соответствующей подобласти имеется главный процессор д0, который не может быть переназначен любой другой подобласти во время вычислительного эксперимента (даже в том случае, когда подобласть не содержит частиц).
Рис. 1. Схема организации параллельных вычислений и балансировки нагрузки между процессорами. Главные процессоры в группе предназначены для вычисления гравитационного потенциала методом свертки. Остальные процессоры в каждой группе занимаются вычислением новых координат частиц.
В результате основной алгоритм выглядит следующим образом:
1. В нулевой момент времени (t = 0):
a. для каждой подобласти вычисляется общее количество частиц, соответствующее числу процессоров в группе и количество частиц, назначенное каждому из процессоров,
b. на каждом процессоре выделяется память для массивов координат и скоростей частиц,
c. выполняется генерация начальных координат и скоростей PIC частиц в области (на каждом процессоре) в соответствии с заданной пользователем функцией распределения.
2. Вычисляется сеточная функция плотности для PIC частиц. Выполняется суммирование плотности на главном процессоре группы: вызывается функция MPI Reduce для процессоров из группы.
3. Вычисляется гравитационный потенциал с помощью параллельного метода свертки.
В том числе выполняются коммуникации МР1_АШоаП между главными процессорами.
4. Для каждого процессора вычисляются новые координаты частиц для следующего шага. Межпроцессорные коммуникации на этом шаге отсутствуют.
5. Вычисляется количество частиц, которое должно находиться в каждой подобласти на следующем шаге. Вычисляется количество процессоров, которое должно быть назначено каждой подобласти.
6. Перераспределяются частицы между процессорными группами. Выполняются межпроцессорные коммуникации МР1_Оа1кегу и МР!^саПету для процессоров из группы.
a. На каждом процессоре могут находиться три типа частиц: (а) частицы, которые остались внутри данной подобласти, (б) частицы, которые должны переместиться в левую смежную подобласть, и (в) частицы, которые должны переместиться в правую смежную подобласть. Общее количество перемещаемых между подобластями частиц невелико из-за необходимости выполнения условия Куранта: оно не может быть больше, чем число частиц, которые содержались в граничных ячейках подобласти. На практике это число существенно меньше и составляет менее 0,1% от общего количества частиц.
b. Вычисляем новое число частиц для каждой из подобластей.
c. Вычисляем число процессоров для каждой подобласти.
^ Определяем те процессоры в каждой группе, которые должны быть переназначены на соседние подобласти (то есть переданы другим группам). Передаем с них те частицы, которые должны остаться в данной подобласти.
е. Передаем частицы, которые перешли из данной подобласти в соседние подобласти.
7. Инкрементируется номер временного шага и выполняется переход на Шаг 2.
МЕТОД СВЕРТКИ ДЛЯ ВЫЧИСЛЕНИЯ ГРАВИТАЦИОННОГО ПОТЕНЦИАЛА
Для вычисления гравитационного потенциала был реализован метод свертки, предложенный Хокни [1]. Его суть заключается в том, что вместо задачи Дирихле для уравнения Пуассона в бесконечной области:
ДФ(г) = 4^Gp(r), Ф(г)||гИм = 0
выполняется эффективное вычисление интеграла, представляющего фундаментальное решение уравнения Пуассона:
Ф(Го)Н^--тт
В дискретном случае в декартовой системе координат на равномерной сетке с числом узлов ^ X Му X и сеточными шагами Ьу, решение будет записываться в следующем виде:
Ф(*о, Уо, 2о) =
=ш
¿=1 _/=1 й = 1
V- хо)2 + О; - Уо)2 + (^л - го)
где = • (^х^у^г) - значения зарядов
(масс), находящихся в узлах сетки.
Воспользовавшись дискретным аналогом теоремы о свертки и алгоритмом быстрого преобразования Фурье, можно вычислить значения Ф(*о,Уо^о) за + +
+ 1о§ операций, что намного быстрее прямого способа вычисления фундаментального решения для уравнения Пуассона и соответствует трудоемкости решения задачи Дирихле для уравнения Пуассона методом разделения переменных [15].
Для того чтобы применить дискретное преобразование Фурье, необходимо устранить неопределенность в точке г = Го и определить функции К и р так, чтобы они стали периодическими.
Это решается с помощью определения сеточной функции К следующим образом:
1
0.5 тт hy,hz)' ВД = 1
{ V х2 +У^2'
|г| = 0, |г| = 0.
Вторая проблема решается с помощью введения фиктивных дополнительных подобластей, дублирующих область решения по каждому направлению в два раза, и доопределения в них функции К так, чтобы она стала периодической во всей новой области, а функция р равной нулю [16].
Параллельная реализация данного алгоритма заключается в разбиении на подобласти и приме-
нении метода транспозиции данных (стандартного подхода для многомерного преобразования Фурье [17]):
1. Вычислительная область в 2D или 3D подразделяется на подобласти в направлении X.
2. Применяется прямое быстрое преобразование Фурье (БПФ) к сеточным функциям плотности ядра потенциала в направлении У и в направлении X.
3. Выполняется транспозиция «слоев» из направления У в направление X.
4. Применяется прямое БПФ в направлении X. Перемножается образ сеточной функции ядра на образ функции плотности. Применяется обратное БПФ в направлении X к полученному результату.
5. Выполняется обратная транспозиция «слоев» из направления X в направление У.
6. Применяется обратное БПФ в направлении X и обратное БПФ в направлении У.
Для выполнения быстрого преобразования Фурье использовалась библиотека FFTW [18].
ТЕСТОВЫЕ РЕЗУЛЬТАТЫ
Тестовые эксперименты проводились на суперкомпьютерах Сибирского суперкомпьютерного центра, Межведомственного суперкомпьютерного центра и суперкомпьютера «Ломоносов» в МГУ. В табл. 1 представлены результаты тестовых экспериментов для параллельного метода свертки в случае плоского двумерного распределения для сеток 16384 X 16384 и 32768 X 32768.
Данный алгоритм обладает хорошей масштабируемостью на большое число процессоров и позволяет вычислять гравитационный потенциал на сетках с числом узлов более 1 млрд за время менее 1 секунды.
Таблица 1
Эксперименты по оценке производительности алгоритма для сеток 16384Х 16384 и 32768 X 32768 при
разном количестве процессоров суперкомпьютера «Ломоносов» МГУ
Число проц. N Время (с) 16384 X16384 Время (с) 32768 X 32768
64 1,75 -
128 1,1 -
256 0,5 2,4
512 0,35 1,3
1024 0,25 0,65
2
В табл. 2 представлены результаты «нагрузочного» эксперимента для расчета одного временного шага. В качестве начального распределения /°(г, и), задавался диск Маклорена с твердотельным вращением и разными дисперсиями скоростей частиц. В качестве технических параметров для полностью трехмерной модели задавалась сетка 1024 X 1024 X 1024, 10 млрд частиц и 1024 процессора. Общее количество подобластей и групп процессоров составляло 128.
Таблица 2 Время расчета одного временного шага для параллельного метода решения системы Власова-Пуассона
Общее время 7 секунд
Particles 60%
MPI Reduce 3%
MPI Alltoall 7,5%
FFT 7,5%
MPI Bcast 7%
MPI_Gatherv, MPI_Scatterv 15%
В табл. 2 относительное разбиение времени по основным программным процедурам и межпроцессорным коммуникациям сгруппировано следующим образом: Particles - вычисление координат частиц (межпроцессорные коммуникации отсутствуют); MPI Reduce - суммирование плотности для подобласти, MPI Alltoall - вычисление гравитационного потенциала; FFT - быстрое преобразование Фурье (межпроцессорные коммуникации отсутствуют); MPI Bcast - пересылка потенциала от главного процессора группы на дочерние; MPI Gatherv, MPI Scatterv - балансировка процессоров и перераспределение частиц.
Табл. 2 показывает, что основная часть времени тратится на последовательные вычисления, в то время как на межпроцессорные коммуникации уходит не более 40%.
ЗАКЛЮЧЕНИЕ
В статье представлен параллельный алгоритм для решения системы уравнений Власова-Пуассона методом частиц. Показано, что данный алгоритм обладает хорошей производительностью и может применяться для суперкомпьютерного моделирования нестационарных задач астрофизики (динамики вращающихся галактик и газопылевых протопланетных дисков), для которых необходимо проведение серийных расчетов с десятками тысяч временных шагов.
СПИСОК ЛИТЕРАТУРЫ
1. Hockney R. W., Eastwood J. W. Computer Simulation Using Particles. New York: McGraw-Hill, 1981. 540 p. [ R. W. Hockney, J. W. Eastwood, Computer Simulation Using Particles. New York: McGraw-Hill, 1981. 540 p. ]
2. Березин Ю. А., Вшивков В. А. Метод частиц в динамике разреженной плазмы. Новосибирск: Наука, 1980. 96 с. [ Yu. A. Berezin and V. A. Vshivkov, Particle-in-cell method in the dynamics of low-density plasma, (in Russian). Novosibirsk: Nauka, 1980. 96 p. ]
3. Barnes J. E., Hut P. A. Hierarchical O(NlogN) Force-Calculation Algorithm // Nature. 1986. Vol. 324, pp. 446-449. [ J.E. Barnes, P.A. Hut, "Hierarchical O(NlogN) Force-Calculation Algorithm", in Nature. 1986. Vol. 324, pp. 446-449. ]
4. Gingold R. A., Monaghan J. J. Smoothed particle hydrodynamics - Theory and application to non-spherical stars // Monthly Notices of the Royal Astronomical Society. 1977. Vol. 181. pp. 375-389. [ R. A. Gingold and J. J. Monaghan, "Smoothed particle hydrodynamics - Theory and application to non-spherical stars", in Monthly Notices of the Royal Astronomical Society. 1977. Vol. 181. pp. 375-389. ]
5. Colella P., Woodward P. R. The Piecewise Parabolic Method (PPM) for gas-dynamical simulations // Journal of Computational Physics. 1984. Vol. 54, No. 1. pp. 174-201. [ P. Colella and P. R. Woodward, "The Piecewise Parabolic Method (PPM) for gas-dynamical simulations", in Journal of Computational Physics. 1984. Vol. 54, No. 1. pp. 174-201. ]
6. Dubeya A., [et al.]. Extensible component-based architecture for FLASH, a massively parallel, multiphysics simulation code // Parallel Computing. 2009. Volume 35. pp. 512-522. [ A. Dubeya, [et al.], "Extensible component-based architecture for FLASH, a massively parallel, multiphysics simulation code", in Parallel Computing. 2009. Volume 35. pp. 512-522. ]
7. Springel V., Yoshida N., White S. D. M. GADGET: a code for collisionless and gasdynamical cosmological simulation // New Astronomy. 2001. Vol. 6. pp. 79-117. [ V. Springel, N. Yo-shida, S. D. M. White, "GADGET: a code for collisionless and gasdynamical cosmological simulation", in New Astronomy. 2001. Vol. 6. pp. 79-117. ]
8. Pearce F. R., Couchman H. M. P. Hydra: a parallel adaptive grid code // New Astronomy. 1997. Vol. 2. No. 5. pp. 411427. [ F. R. Pearce and H. M. P. Couchman, "Hydra: a parallel adaptive grid code", in New Astronomy. 1997. Vol. 2. No. 5. pp. 411-427. ]
9. Springel V. [et al]. Simulations of the formation, evolution and clustering of galaxies and quasars // Nature. 2005. Vol. 435. pp. 629-636. [ V. Springel, [et al.], "Simulations of the formation, evolution and clustering of galaxies and quasars", in Nature. 2005. Vol. 435. pp. 629-636. ]
10. Klypin A. A. [et al.]. Dark Matter Halos in the Standard Cosmological Model: Results from the Bolshoi Simulation // As-trophysical Journal. 2011. Vol. 740. pp. 102-118. [ A. A. Klypin, [et al.], "Dark Matter Halos in the Standard Cosmological Model: Results from the Bolshoi Simulation", in Astrophysical Journal. 2011. Vol. 740. pp. 102-118. ]
11. Feng Y., [et al.]. BlueTides: First galaxies and reioniza-tion // Monthly Notice of Royal Astronomical Society. 2015 (submitted). [ Y. Feng, [et al], "BlueTides: First galaxies and reionization", in Monthly Notice of Royal Astronomical Society. 2016 (submitted) ]
12. Снытников Н. В., Вшивков В. А. Метод декомпозиции области для суперкомпьютерного моделирования гра-витирующих систем // Суперкомпьютерные дни в России: международная конференция (Москва, 28-29 сентября, 2015): тр. конф. 2015. С. 572-580. [ N. V. Snytnikov and V. A. Vshivkov, "Domain decomposition method for superco-muter simulation of gravitating systems", (in Russian), in Proc. Int. Conference "Supercomputing Days in Russia", Moscow, Russia, 2015, pp.572-580. ]
13. Снытников В. Н. [и др.]. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астрономический журнал. 2004. Т. 30, №.2. С. 146-160 [V. N. Snytnikov, [et al.] "Three-Dimensional Numerical Simulation of a Nonstationary Gravitating N-Body System with Gas", (in Russian), in Pisma v Astromo-michieskiiZhurnal. 2004. Vol. 30. P. 124-137. ]
14. Вшивков В. А., Снытников В. Н., Снытников Н. В. Моделирование трехмерной динамики вещества в гравитационном поле на многопроцессорных ЭВМ // Вычислительные технологии. 2006, Т.11. № 2. С. 15-27. [ V. A. Vshivkov, V. N. Snytnikov, N. V. Snytnikov. "Simulation of three-dimensional dynamics of matter in gravitational field with the use of multiprocessor computer", (in Russian), in Vychislitelnye tekhnologii. 2006. Vol. 11, No. 2. pp. 15-27. ]
15. Самарский А. А., Андреев В. Б. Разностные методы для эллиптических уравнений. М.: Наука, 1976. 352 c. [ A. A. Samarskii, V. B. Andreev, Difference methods for elliptic equations, (in Russian). Moscow: Nauka, 1976. p. 352. ]
16. Eastwood J. W., Brownrigg D. R. K. Remarks on the Solution of Poisson's Equation for Isolated Systems // Journal of Computational Physics. 1979. Vol. 32. pp. 24-38. [ J. W. Eastwood and D. R. K. Brownrigg, "Remarks on the Solution of Pois-son's Equation for Isolated Systems", in Journal of Computational Physics. 1979. Vol. 32. pp. 24-38. ]
17. Ayala O., Wang L. P. Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition // Parallel Computing. 2013. Vol.39. pp. 58-77. [ O. Ayala and L. P. Wang, "Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition", in Parallel Computing. 2013. Vol.39. pp. 58-77. ]
18. Frigo M., Johnson S.G. FFTW software. [Электронный ресурс]. URL: http://www.fftw.org (дата обращения: 01.02.2016). [ M. Frigo, and S. G. Johnson (2016, Feb. 12), "FFTW software" [Online]. Available: http://www.fftw.org ]
ОБ АВТОРЕ
СНЫТНИКОВ Николай Валерьевич, науч. сотр. лаборатории параллельных алгоритмов решения больших задач Института вычислительной математики и математической геофизики СО РАН. Иссл. в обл. суперкомпьютерного моделирования астрофизических систем.
METADATA
Title: Scalable parallel algorithm for simulation of 3D dynamics
of gravitating systems using particle-in-cells methods. Author: N.V. Snytnikov12 Affiliation:
1 Institute of Computational Mathematics and Mathematical Geophysics, Russia.
2 Novosibirsk State University, Russia Email: 1 [email protected] Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 20, no. 3 (73), pp. 137-142, 2016. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: Parallel algorithm for numerical solving of Vlasov-Poisson equation using particle-in-cell method is developed. New method of processors dynamic load-balancing is used. Processors are distributed among subdomains corresponding to the number of particles. Domain decomposition method allows to combine grid (Eulerian) method for solving Poisson equation and Lagrangian method for solving Vlasov equation. It uses physical features of problems of non-stationary rotating gravitating disks (2D and 3D). Key words: gravitating systems, Poisson equation, Vlasov equation, particle-in-cell method, dynamic load balancing About author:
SNYTNIKOV Nikolay Valerievich, scientific researcher in the Laboratory of Parallel Algorithms for Solving Big Problems in Institute of Computational Mathematics and Mathematical Geophysics. Researcher in the field of supercomputer simulation of astrophysical systems.
