Предложен параллельный алгоритм для численного решения системы уравнений Власова-Пуассона методом частиц. Используется новый метод динамической балансировки процессоров, которые распределяются между подобластями в соответствии с реальным числом модельных частиц. Метод декомпозиции области позволяет комбинировать сеточный (эйлеров) метод для решения уравнения Пуассона с лагранжевым методом частиц для решения уравнения Власова. Он учитывает физические особенности задач моделирования нестационарных вращающихся дисков (двумерных и трехмерных). Ключевые слова: гравитирующие системы, уравнение Пуассона, уравнение Власова, метод частиц, динамическая балансировка.
и уравнения Пуассона для гравитационного потенциала \Phi = \Phi( t, \bfr ) с краевыми условиями
D\eltaPhi( t, \bfr ) = 4\piG\rho (t, \bfr ),
P\hi( t, \bfr )
= 0, | | \bfr грант № 14-11-00485. Программная реализация алгоритма и тестовые эксперименты выполнены при подagora.guru.ru/pavt где G – гравитационная постоянная. ее с функцией распределения f : Замыкает систему следующее уравнение для плотности \rho = \rho (t, \bfr ), которое связывает i\nt
Для решения уравнения Власова используется метод частиц в ячейках [
1024 \times В типичном случае максимальное количество частиц, которое можно поместить в одно вычислительное устройство (CPU, GPU, Intel Phi), составляет 8-16 миллионов, а максимальное количество узлов сетки составляет 256 \times 256 \times 256 для трехмерных задач и 4096 для двумерных. Это означает, что большие вычислительные области (такие 1024) должны подразделяться на меньшие подобласти. В то же самое время модельные частицы могут многократно перелетать из одной области в другую в течение одного численного эксперимента (например, если рассматривается моделирование дисков на временных масштабах порядка десятков оборотов, где каждая частица двигается по эллиптической или эпициклической орбите вокруг общего центра масс), что создает проблему пересылок данных на каждом временном шаге. Кроме того, частицы могут быть распределены неравномерно между подобластями. Например, во вращающемся диске могут возникать гравитационные неустойчивости, которые приводят к появлению разнообразных структур — сгустков вещества или спиральных волн, которые способны передвигаться по всей вычислительной подобласти. Плотность вещества (и количество частиц) в этих структурах может быть на порядки больше фоновых значений. Следовательно, количество частиц в каждой подобласти может так же отличаться на порядки и меняться со временем, особенно в тех случаях, когда сгусток, содержащий большое количество частиц, двигается Разработанный параллельный алгоритм использует следующие свойства задачи и чиссетки в типичных расчетах находится в диапазоне 10 \div 100.
подобластями).
b\ulet Чтобы удовлетворять условию Куранта устойчивости численного метода, частица
может перелететь за один шаг не далее, чем в соседнюю ячейку. Это означает, что
число частиц, которые должны перемещаться между смежными подобластями, намного
меньше, чем общее число частиц в смежных подобластях (порядка 10 \cdot ny для 2D
задач и 10 \cdot nynz для 3D задач, где ny, nz — это число всех граничных узлов между
Общая схема параллельного алгоритма описана в разделе 2, а детальное описание
алгоритмов приведено в статье [
Исходная вычислительная область (2D или 3D) равномерно подразделяется на K
подобластей одинакового размера S1, S2, .., SK в направлении X. Группа процессоров Gk
назначается каждой подобласти Sk (см. Рис. 1). Количество процессоров в группе Gk равно
Pk, Pk \geq 1. Общее количество процессоров равно P = \sum Pk. Количество узлов (ячеек) в
каждой вычислительной подобласти выбирается таким образом, что все сеточные функции
(гравитационный потенциал, плотность, сила) должны полностью помещаться в
оперативную память одного вычислительного устройства и оставить в памяти достаточно места для
хранения массивов с частицами (т.е. их пространственных координат и скоростей). Обычно
это число равно 0.5 \div 2 миллионов узлов (соответствует сетке 1283 или 10242). Количество
процессоров в группе Gk выбирается так, чтобы обеспечить приблизительно одинаковое
количество частиц на одном процессоре (см. Алгоритм 2 в [
В результате основной алгоритм выглядит следующим образом: 1. В нулевой момент времени (t = 0): (a) для кажой подобласти Sk вычисляется общее количество частиц, соответствующее числу процессоров в группе Gk, и количество частиц, назначенное каждому из процессоров, (b) на каждом процессоре выделяется память для массивов координат и скоростей частиц, (c) выполняется генерация начальных координат и скоростей PIC частиц в области (на каждом процессоре) в соответствии с заданной пользователем функцией распределения. agora.guru.ru/pavt ние плотности на главном процессоре группы: вызывается функция MPI_Reduce для процессоров из группы. том числе выполняются коммуникации MPI_Alltoall между главными процессорами. 4. Для каждого процессора вычисляются новые координаты частиц для следующего шага. Межпроцессорные коммуникации на этом шаге отсутствуют. 6. Перераспределяются частицы между процессорными группами. Выполняются межпроцессорные коммуникации MPI_Gatherv и MPI_Scatterv для процессоров из группы. (a) На каждом процессоре могут находиться три типа частиц: (а) частицы, которые остались внутри данной подобласти, (б) частицы, которые должны переместиться в левую смежную подобласть, и (в) частицы, которые должны переместиться в правую смежную подобласть. Общее количество перемещаемых между подобластями частиц невелико из-за необходимости выполнения условия Куранта: оно не может быть больше, чем число частиц, которые содержались в граничных ячейках подобласти. На практике это число существенно меньше и составляет менее 0.1% от общего количества частиц. (b) Вычисляем новое число частиц для каждой из подобластей. (c) Вычисляем число процессоров для каждой подобласти. (d) Определяем те процессоры в каждой группе, которые должны быть переназначены на соседние подобласти (то есть переданы другим группам). Передаем с них те частицы, которые должны остаться в данной подобласти. (e) Передаем частицы, которые перешли из данной подобласти в соседние подобласти.
7. Инкрементируется номер временного шага и выполняется переход на Шаг 2. 3. Метод свёртки для вычисления гравитационного потенциала Для вычисления гравитационного потенциала был реализован метод свёртки, предложенный Хокни [1]. Его суть заключается в том, что вместо задачи Дирихле для уравнения Пуассона в бесконечной области: выполняется эффективное вычисление интеграла, представляющего фундаментальное решение уравнения Пуассона: В дискретном случае в декартовой системе координат на равномерной сетке c числом узлов Nx \times
Ny \times Nz и сеточными шагами hx, hy, hz решение будет записываться в следующем виде: P\hi( x0, y0, z0) = s\um Nx- 1 Ny- 1 Nz- 1 s\um
s\um s\qrt{} (xi - x0)2 + (yi - y0)2 + (zi - z0)2 agora.guru.ru/pavt где qi,j,k = \rho i,j,k \cdot (hxhyhz) — значения зарядов (масс), находящихся в узлах сетки.
Воспользовавшись дискретным аналогом теоремы о свёртки и алгоритмом быстрого
преобразования Фурье, можно вычислить сумму (3) за O(NxNyNz(log2Nx+log2Ny+log2Nz))
операций, что намного быстрее прямого способа вычисления фундаментального решения
для уравнения Пуассона и соответствует трудоемкости решения задачи Дирихле для
уравнения Пуассона методом разделения переменных [
Для того, чтобы применить дискретное преобразование Фурье, необходимо устранить неопределенность в точке \bfx
= \bfx \bfzero и определить функции K и \rho так, чтобы они стали периодическими. Это решается с помощью определения сеточной функции K следующим образом:
K(x, y, z) = l\eft{ 0.5 \mathr{}in(
1 s\urd x2+y2+z2 hx,hy,hz) , \sqrt{} x2 + y2 + z2 = 0, 1
, \sqrt{} x2 + y2 + z2 > 0.
Вторая проблема решается с помощью введения фиктивных дополнительных
подобластей, дублирующих область решения по каждому направлению в два раза, и доопределения
в них функции K так, чтобы она стала периодической во всей новой области, а функция
r\ho — равной нулю [
Параллельная реализация данного алгоритма заключается в разбиении на подобласти
и применении метода транспозиции данных (стандартного подхода для многомерного
преобразования Фурье [
Для выполнения быстрого преобразования Фурье использовалась библиотека FFTW [
Тестовые эксперименты проводились на суперкомпьютерах Сибирского суперкомпьютерного центра, Межведомственного суперкомпьютерного центра и суперкомпьютера «Ломоносов» в МГУ. В Таблице 1 представлены результаты тестовых экспериментов для параллельного метода свертки в случае плоского двумерного распределения для сеток 16384\times 16384 и 32768\times 32768. Данный алгоритм обладает хорошей масштабируемостью на большое число процессоров и позволяет вычислять гравитационный потенциал на сетках с числом узлов более 1 млрд. за время менее 1 секунды.
В Таблице 2 представлены результаты «нагрузочного» эксперимента для расчета одного временного шага. В качестве начального распределения f 0(\bfr , \bfu ) задавался диск Маклорена с твердотельным вращением и разными дисперсиями скоростей частиц. В качестве технических параметров для полностью трехмерной модели задавалась сетка 1024 \times 1024 \times 1024, 0.65 Таблица 1. Эксперименты по оценке производительности алгоритма для сетки 16384 \times 16384 и 32768 \times 32768 при разном количестве процессоров суперкомпьютера «Ломоносов» МГУ. Число процессоров N Время (c) 16384 \times 16384
Время (c) 32768 \times 32768 Таблица 2. Время расчета одного временного шага для параллельного метода решения системы Власова-Пуассона. Относительное разбиение времени по основным программным процедурам и межпроцессорным коммуникациям сгруппировано следующим образом: Particles: вычисление координат частиц (межпроцессорные коммуникации отсутствуют), MPI_Reduce: суммирование плотности для подобласти, MPI_Alltoall: вычисление гравитационного потенциала, FFT: быстрое преобразование Фурье (межпроцессорные коммуникации отсутствуют), MPI_Bcast: пересылка потенциала от главного процессора группы на дочерние, MPI_Gatherv, MPI_Scatterv: балансировка процессоров и перераспределение частиц.
Общее время
MPI_Gatherv, MPI_Scatterv 7 секунд 60% 3% 7.5% 7.5% 7% 15% 10 миллиардов частиц и 1024 процессора. Общее количество подобластей и групп процессоров составляло 128.
Таблица показывает, что основная часть времени тратится на последовательные вычисления, в то время как на межпроцессорные коммуникации уходит не более 40%. 5. Заключение
В статье представлен параллельный алгоритм для решения системы уравнений ВласоваПуассона методом частиц. Показано, что данный алгоритм обладает хорошей производительностью и может применяться для суперкомпьютерного моделирования нестационарных задач астрофизики (динамики вращающихся галактик и газопылевых протопланетных дисков), для которых необходимо проведение серийных расчетов с десятками тысяч временных шагов. Литература 1. Hockney, R.W., Eastwood, J.W. Computer Simulation Using Particles. McGraw-Hill. New
York. 1981. 540 p. 3. Barnes J.E., Hut P. A Hierarchical O(NlogN) Force-Calculation Algorithm // Nature.
1986. Vol. 324, P. 446-449. 4. Gingold R.A. and Monaghan J.J. Smoothed particle hydrodynamics - Theory and application to non-spherical stars // Monthly Notices of the Royal Astronomical Society. 1977. Vol. 181. P. 375-389. 5. Colella P., Woodward P.R. The Piecewise Parabolic Method (PPM) for gas-dynamical simulations // Journal of Computational Physics. 1984. Vol. 54, No. 1. P. 174-201. 6. Dubeya A., Antypasb K., Ganapathyc M.K., et al. Extensible component-based architecture for FLASH, a massively parallel, multiphysics simulation code // Parallel Computing. 2009. Volume 35. P. 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. P. 79-117. 8. Pearce F.R., Couchman H.M.P. Hydra: a parallel adaptive grid code // New Astronomy.
1997. Vol. 2, No. 5. P. 411-427. 9. Springel V. et al. Simulations of the formation, evolution and clustering of galaxies and quasars // Nature. 2005. Vol. 435. P. 629. 10. Klypin A.A. et al. Dark Matter Halos in the Standard Cosmological Model: Results from the Bolshoi Simulation // Astrophysical Journal. 2011. Vol. 740. P. 102. 11. Feng Y., Di Matteo T., Croft R.A.C., Bird S., Battaglia N., Wilkins S. BlueTides: First galaxies and reionization // Monthly Notice of Royal Astronomical Society. 2015 (submitted) 13. Снытников В.Н., Вшивков В.А., Кукшева Э.А., Неупокоев Е.В., Никитин С.А., Снытников А.В. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астрономический журнал. 2004. Т. 30, №.2. 146-160.
Scalable parallel algorithm for simulation of 3D dynamics of gravitating systems using particle-in-cell method
N.V. Snytnikov1 Institute of Computational Mathematics and Mathematical Geophysics 1 We developed parallel algorithm for solving Vlasov-Poisson systems of equations using particle-in-cell method. It uses new technique of dynamic load balancing for processors, which are distributed between subdomains in correspndance with the number of modeling particles located in the subdomain. Domain decomposition method combines grid (eulerian) method for solving Poisson equation with lagrangian paricle method for solving Vlasov equation. It takes into account physical features of the modeling nonstationary rotating disks (both 2D and 3D). 1. Hockney, R.W., Eastwood, J.W. Computer Simulation Using Particles. McGraw-Hill. New
York. 1981. 540 p. 3. Barnes J.E., Hut P. A Hierarchical O(NlogN) Force-Calculation Algorithm // Nature.
1986. Vol. 324, P. 446-449. 4. Gingold R.A. and Monaghan J.J. Smoothed particle hydrodynamics - Theory and application to non-spherical stars // Monthly Notices of the Royal Astronomical Society. 1977. Vol. 181. P. 375-389. 5. Colella P., Woodward P.R. The Piecewise Parabolic Method (PPM) for gas-dynamical simulations // Journal of Computational Physics. 1984. Vol. 54, No. 1. P. 174-201. 6. Dubeya A., Antypasb K., Ganapathyc M.K., et al. Extensible component-based architecture for FLASH, a massively parallel, multiphysics simulation code // Parallel Computing. 2009. Volume 35. P. 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. P. 79-117. 8. Pearce F.R., Couchman H.M.P. Hydra: a parallel adaptive grid code // New Astronomy.
1997. Vol. 2, No. 5. P. 411-427. 9. Springel V. et al. Simulations of the formation, evolution and clustering of galaxies and quasars // Nature. 2005. Vol. 435. P. 629. 10. Klypin A.A. et al. Dark Matter Halos in the Standard Cosmological Model: Results from the Bolshoi Simulation // Astrophysical Journal. 2011. Vol. 740. P. 102. 11. Feng Y., Di Matteo T., Croft R.A.C., Bird S., Battaglia N., Wilkins S. BlueTides: First galaxies and reionization // Monthly Notice of Royal Astronomical Society. 2015 (submitted)