Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Масштабируемый параллельный алгоритм для моделирования трехмерной динамики \ast гравитирующих систем методом частиц 1 Н.В. Снытников 1 Институт Вычислительной Математики и Математической Геофизики СО РАН Предложен параллельный алгоритм для численного решения системы уравне- ний Власова-Пуассона методом частиц. Используется новый метод динамиче- ской балансировки процессоров, которые распределяются между подобластями в соответствии с реальным числом модельных частиц. Метод декомпозиции обла- сти позволяет комбинировать сеточный (эйлеров) метод для решения уравнения Пуассона с лагранжевым методом частиц для решения уравнения Власова. Он учитывает физические особенности задач моделирования нестационарных вра- щающихся дисков (двумерных и трехмерных). Ключевые слова: гравитирующие системы, уравнение Пуассона, уравнение Вла- сова, метод частиц, динамическая балансировка. 1. Введение Суперкомпьютерное моделирование динамики гравитирующих систем (таких как га- лактики или околозвездные диски) является одной из наиболее трудных и актуальных проблем в современной вычислительной астрофизике. Хотя однопроцессорные численные методы были созданы и успешно используются уже более 30 лет [1–5], разработка соответ- ствующих им параллельных алгоритмов требует значительных модификаций, связанных с адаптацией к новым аппаратным и програмнным архитектурам суперкомпьютеров. Другие трудности возникают из-за разнообразия моделируемых астрофизических процессов [6–8], для которых существуют специализированные программные пакеты [9–11] не всегда при- менимые для решения смежных задач. В данной статье предложен параллельный алгоритм для моделирования вращающихся бесстолкновительных гравитирующих систем методом частиц в ячейках, и приведены ре- зультаты численных экспериментов (детальное описание алгоритмов приведено в [12]). Ал- горитм может быть применен как для модели тонкого диска (когда отсутствует движение вещества в вертикальном направлении), так и для полностью трехмерной модели [13, 14]. Новым в этом подходе является метод динамической балансировки загрузки процессоров при вычислении траектории модельных частиц, которые в случае вращающихся систем мо- гут многократно переходить между подобластями (и соответствующим им процессорам). Решаемая система уравнений состоит из уравнения Власова для функции распреде- ления вещества f = f (t, \bfr , \bfu ), зависящей от времени t, пространственной координаты \bfr и вектора скорости \bfu , с заданным начальным значением f 0 (\bfr , \bfu ) (в виде некоторого враща- ющегося диска): \partial f \partial f \partial f + \bfu - \nabla \Phi (t, \bfr ) = 0, f (0, \bfr , \bfu ) = f 0 (\bfr , \bfu ) \partial t \partial \bfr \partial \bfu и уравнения Пуассона для гравитационного потенциала \Phi = \Phi (t, \bfr ) с краевыми условиями для изолированных систем: \Delta \Phi (t, \bfr ) = 4\pi G\rho (t, \bfr ), \Phi (t, \bfr )| | \bfr | \rightarrow \infty = 0, \ast Разработка теоретических основ алгоритма выполнена при поддержке Российского Научного Фонда, грант № 14-11-00485. Программная реализация алгоритма и тестовые эксперименты выполнены при под- держке РФФИ, гранты № 16-07-00916, № 14-07-00241, № 14-01-00392. 298 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt где G – гравитационная постоянная. Замыкает систему следующее уравнение для плотности \rho = \rho (t, \bfr ), которое связывает ее с функцией распределения f : \int \rho (t, \bfr ) = f (t, \bfr , \bfu )d\bfu . \bfu Для решения уравнения Власова используется метод частиц в ячейках [1, 2], а для ре- шения уравнения Пуассона — метод свёртки. Соответствующий параллельный алгоритм должен подразделять исходную вычислительную область на такие подобласти, что каждая из них могла бы быть обработана одним вычислительным устройством с 1-4 ГБ памяти, и распределять модельные частицы между процессорами так, чтобы они могли иметь до- ступ к значениям требуемых сеточных функций в подобластях, свободно двигаться между подобластями, и общее количество частиц, передаваемых между процессорами, было ми- нимальным. В типичном случае максимальное количество частиц, которое можно поместить в од- но вычислительное устройство (CPU, GPU, Intel Phi), составляет 8-16 миллионов, а мак- симальное количество узлов сетки составляет 256 \times 256 \times 256 для трехмерных задач и 4096 \times 4096 для двумерных. Это означает, что большие вычислительные области (такие как 1024 \times 1024 \times 1024) должны подразделяться на меньшие подобласти. В то же самое время модельные частицы могут многократно перелетать из одной области в другую в течение одного численного эксперимента (например, если рассматривается моделирование дисков на временных масштабах порядка десятков оборотов, где каждая частица двигает- ся по эллиптической или эпициклической орбите вокруг общего центра масс), что создает проблему пересылок данных на каждом временном шаге. Кроме того, частицы могут быть распределены неравномерно между подобластями. Например, во вращающемся диске могут возникать гравитационные неустойчивости, которые приводят к появлению разнообразных структур — сгустков вещества или спиральных волн, которые способны передвигаться по всей вычислительной подобласти. Плотность вещества (и количество частиц) в этих струк- турах может быть на порядки больше фоновых значений. Следовательно, количество ча- стиц в каждой подобласти может так же отличаться на порядки и меняться со временем, особенно в тех случаях, когда сгусток, содержащий большое количество частиц, двигается вокруг центра масс. Разработанный параллельный алгоритм использует следующие свойства задачи и чис- ленного метода: \bullet Трудоемкость решения уравнения Власова (т.е. интегрирование траекторий PIC-частиц) существенно выше (от 10 до 100 раз), чем трудоемкость вычисления гравитационного потенциала. Это связано с тем, что среднее число модельных частиц на одну ячейку сетки в типичных расчетах находится в диапазоне 10 \div 100. \bullet Чтобы удовлетворять условию Куранта устойчивости численного метода, частица мо- жет перелететь за один шаг не далее, чем в соседнюю ячейку. Это означает, что чис- ло частиц, которые должны перемещаться между смежными подобластями, намного меньше, чем общее число частиц в смежных подобластях (порядка 10 \cdot ny для 2D за- дач и 10 \cdot ny nz для 3D задач, где ny , nz — это число всех граничных узлов между подобластями). Общая схема параллельного алгоритма описана в разделе 2, а детальное описание ал- горитмов приведено в статье [12]. Раздел 3 посвящен описанию метода вычисления грави- тационного потенциала. В разделе 4 приводятся результаты экспериментов. 299 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 2. Параллельный алгоритм Исходная вычислительная область (2D или 3D) равномерно подразделяется на K под- областей одинакового размера S1 , S2 , .., SK в направлении X. Группа процессоров Gk на- значается каждой подобласти Sk (см. Рис. 1). Количество процессоров в группе Gk равно \sum Pk , Pk \geq 1. Общее количество процессоров равно P = Pk . Количество узлов (ячеек) в каждой вычислительной подобласти выбирается таким образом, что все сеточные функции (гравитационный потенциал, плотность, сила) должны полностью помещаться в оператив- ную память одного вычислительного устройства и оставить в памяти достаточно места для хранения массивов с частицами (т.е. их пространственных координат и скоростей). Обычно это число равно 0.5 \div 2 миллионов узлов (соответствует сетке 1283 или 10242 ). Количество процессоров в группе Gk выбирается так, чтобы обеспечить приблизительно одинаковое количество частиц на одном процессоре (см. Алгоритм 2 в [12]). Это означает, что если под- область Sk1 содержит больше частиц, чем подобласть Sk2 , то число процессоров Pk1 в группе Gk1 будет больше или равно количеству Pk2 в группе Gk2 . В каждой группе процессоров Gk и соответствующей подобласти Sk имеется главный процессор gk0 , который не может быть переназначен любой другой подобласти во время вычислительного эксперимента (даже в том случае, когда подобласть Sk не содержит частиц). Рис. 1. Схема организации параллельных вычислений и балансировки нагрузки между процессо- рами. Главные процессоры в группе предназначены для вычисления гравитационного потенциала методом свёртки. Остальные процессоры в каждой группе занимаются вычислением новых коор- динат частиц В результате основной алгоритм выглядит следующим образом: 1. В нулевой момент времени (t = 0): (a) для кажой подобласти Sk вычисляется общее количество частиц, соответствую- щее числу процессоров в группе Gk , и количество частиц, назначенное каждому из процессоров, (b) на каждом процессоре выделяется память для массивов координат и скоростей частиц, (c) выполняется генерация начальных координат и скоростей PIC частиц в обла- сти (на каждом процессоре) в соответствии с заданной пользователем функцией распределения. 300 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 2. Вычисляется сеточная функция плотности для PIC частиц. Выполняется суммирова- ние плотности на главном процессоре группы: вызывается функция MPI_Reduce для процессоров из группы. 3. Вычисляется гравитационный потенциал с помощью параллельного метода свертки. В том числе выполняются коммуникации MPI_Alltoall между главными процессорами. 4. Для каждого процессора вычисляются новые координаты частиц для следующего ша- га. Межпроцессорные коммуникации на этом шаге отсутствуют. 5. Вычисляется количество частиц, которое должно находиться в каждой подобласти Sk на следующем шаге. Вычисляется количество процессоров, которое должно быть назначено каждой подобласти. 6. Перераспределяются частицы между процессорными группами. Выполняются меж- процессорные коммуникации MPI_Gatherv и MPI_Scatterv для процессоров из груп- пы. (a) На каждом процессоре могут находиться три типа частиц: (а) частицы, которые остались внутри данной подобласти, (б) частицы, которые должны переместить- ся в левую смежную подобласть, и (в) частицы, которые должны переместиться в правую смежную подобласть. Общее количество перемещаемых между подоб- ластями частиц невелико из-за необходимости выполнения условия Куранта: оно не может быть больше, чем число частиц, которые содержались в граничных ячейках подобласти. На практике это число существенно меньше и составляет менее 0.1% от общего количества частиц. (b) Вычисляем новое число частиц для каждой из подобластей. (c) Вычисляем число процессоров для каждой подобласти. (d) Определяем те процессоры в каждой группе, которые должны быть переназначе- ны на соседние подобласти (то есть переданы другим группам). Передаем с них те частицы, которые должны остаться в данной подобласти. (e) Передаем частицы, которые перешли из данной подобласти в соседние подобла- сти. 7. Инкрементируется номер временного шага и выполняется переход на Шаг 2. 3. Метод свёртки для вычисления гравитационного потенциала Для вычисления гравитационного потенциала был реализован метод свёртки, предло- женный Хокни [1]. Его суть заключается в том, что вместо задачи Дирихле для уравнения Пуассона в бесконечной области: \Delta \Phi (\bfx ) = \rho (\bfx ), (1) \Phi (\bfx )| \bfx \rightarrow \infty = 0, выполняется эффективное вычисление интеграла, представляющего фундаментальное ре- шение уравнения Пуассона: \int \rho (\bfx )d\bfx \Phi (\bfx \bfzero ) = - . (2) | \bfx \bfzero - \bfx | В дискретном случае в декартовой системе координат на равномерной сетке c числом узлов Nx \times Ny \times Nz и сеточными шагами hx , hy , hz решение будет записываться в следующем виде: y - 1 Nz - 1 x - 1 N\sum N\sum \sum qi,j,k \Phi (x0 , y0 , z0 ) = - \sqrt{} , (3) i=1 j=1 k=1 (xi - x0 ) + (yi - y0 )2 + (zi - z0 )2 2 301 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt где qi,j,k = \rho i,j,k \cdot (hx hy hz ) — значения зарядов (масс), находящихся в узлах сетки. Воспользовавшись дискретным аналогом теоремы о свёртки и алгоритмом быстрого преобразования Фурье, можно вычислить сумму (3) за O(Nx Ny Nz (log2 Nx +log2 Ny +log2 Nz )) операций, что намного быстрее прямого способа вычисления фундаментального решения для уравнения Пуассона и соответствует трудоемкости решения задачи Дирихле для урав- нения Пуассона методом разделения переменных [15]. Для того, чтобы применить дискретное преобразование Фурье, необходимо устранить неопределенность в точке \bfx = \bfx \bfzero и определить функции K и \rho так, чтобы они стали пе- риодическими. Это решается с помощью определения сеточной функции K следующим образом: \left\{ 1 \sqrt{} 0.5 \mathrm{m}\mathrm{i}\mathrm{n}(hx,hy,hz) , x2 + y 2 + z 2 = 0, K(x, y, z) = \sqrt{} \surd 1 , x2 + y 2 + z 2 > 0. x2 +y 2 +z 2 Вторая проблема решается с помощью введения фиктивных дополнительных подобла- стей, дублирующих область решения по каждому направлению в два раза, и доопределения в них функции K так, чтобы она стала периодической во всей новой области, а функция \rho — равной нулю [16]. Параллельная реализация данного алгоритма заключается в разбиении на подобласти и применении метода транспозиции данных (стандартного подхода для многомерного пре- образования Фурье [18]): 1. Вычислительная область в 2D или 3D подразделяется на подобласти в направлении X. 2. Применяется прямое быстрое преобразование Фурье (БПФ) к сеточным функциям плотности ядра потенциала в направлении Y и в направлении Z. 3. Выполняется транспозиция «слоев» из направления Y в направление X. 4. Применяется прямое БПФ в направлении X. Перемножается образ сеточной функции ядра на образ функции плотности. Применяется обратное БПФ в направлении X к полученному результату. 5. Выполняется обратная транспозиция «слоев» из направления X в направление Y. 6. Применяется обратное БПФ в направлении Z и обратное БПФ в направлении Y. Для выполнения быстрого преобразования Фурье использовалась библиотека FFTW [17]. 4. Тестовые результаты Тестовые эксперименты проводились на суперкомпьютерах Сибирского суперкомпью- терного центра, Межведомственного суперкомпьютерного центра и суперкомпьютера «Ло- моносов» в МГУ. В Таблице 1 представлены результаты тестовых экспериментов для па- раллельного метода свертки в случае плоского двумерного распределения для сеток 16384\times 16384 и 32768\times 32768. Данный алгоритм обладает хорошей масштабируемостью на большое число процессоров и позволяет вычислять гравитационный потенциал на сетках с числом узлов более 1 млрд. за время менее 1 секунды. В Таблице 2 представлены результаты «нагрузочного» эксперимента для расчета одного временного шага. В качестве начального распределения f 0 (\bfr , \bfu ) задавался диск Маклорена с твердотельным вращением и разными дисперсиями скоростей частиц. В качестве техни- ческих параметров для полностью трехмерной модели задавалась сетка 1024 \times 1024 \times 1024, 302 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Таблица 1. Эксперименты по оценке производительности алгоритма для сетки 16384 \times 16384 и 32768 \times 32768 при разном количестве процессоров суперкомпьютера «Ломоносов» МГУ. Число процессоров N Время (c) 16384 \times 16384 Время (c) 32768 \times 32768 64 1.75 — 128 1.1 — 256 0.5 2.4 512 0.35 1.3 1024 0.25 0.65 Таблица 2. Время расчета одного временного шага для параллельного метода решения систе- мы Власова-Пуассона. Относительное разбиение времени по основным программным процедурам и межпроцессорным коммуникациям сгруппировано следующим образом: Particles: вычисление коор- динат частиц (межпроцессорные коммуникации отсутствуют), MPI_Reduce: суммирование плотно- сти для подобласти, MPI_Alltoall: вычисление гравитационного потенциала, FFT: быстрое преобра- зование Фурье (межпроцессорные коммуникации отсутствуют), MPI_Bcast: пересылка потенциала от главного процессора группы на дочерние, MPI_Gatherv, MPI_Scatterv: балансировка процессо- ров и перераспределение частиц. Общее время 7 секунд Particles 60% MPI_Reduce 3% MPI_Alltoall 7.5% FFT 7.5% MPI_Bcast 7% MPI_Gatherv, MPI_Scatterv 15% 303 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 10 миллиардов частиц и 1024 процессора. Общее количество подобластей и групп процес- соров составляло 128. Таблица показывает, что основная часть времени тратится на последовательные вычис- ления, в то время как на межпроцессорные коммуникации уходит не более 40%. 5. Заключение В статье представлен параллельный алгоритм для решения системы уравнений Власова- Пуассона методом частиц. Показано, что данный алгоритм обладает хорошей производи- тельностью и может применяться для суперкомпьютерного моделирования нестационарных задач астрофизики (динамики вращающихся галактик и газопылевых протопланетных дис- ков), для которых необходимо проведение серийных расчетов с десятками тысяч временных шагов. Литература 1. Hockney, R.W., Eastwood, J.W. Computer Simulation Using Particles. McGraw-Hill. New York. 1981. 540 p. 2. Березин Ю.А., Вшивков В.А. Метод частиц в динамике разреженной плазмы. Наука. Новосибирск. 1980. 96 с. 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) 12. Снытников Н.В., Вшивков В.А. Метод декомпозиции области для суперкомпьютерного моделирования гравитирующих систем // СУПЕРКОМПЬЮТЕРНЫЕ ДНИ В РОССИИ. Труды международной конференции. Под редакцией чл.-корр. РАН Вл.В. Воеводина. 2015. С. 572-580. 304 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 13. Снытников В.Н., Вшивков В.А., Кукшева Э.А., Неупокоев Е.В., Никитин С.А., Снытников А.В. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астрономический журнал. 2004. Т. 30, №.2. 146-160. 14. Вшивков В.А., Снытников В.Н., Снытников Н.В. Моделирование трехмерной динамики вещества в гравитационном поле на многопроцессорных ЭВМ // Вычислительные технологии. 2006, Т.11. N.2. С.15-27. 15. Самарский А.А., Андреев В.Б. Разностные методы для эллиптических уравнений. М.: Наука, 1976. 352 c. 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, P.24-38. 17. M. Frigo, and S.G. Johnson. FFTW software. URL: http://www.fftw.org (дата обращения: 01.02.2016) 18. Ayala O., Wang L.P. Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition // Parallel Computing. 2013. Vol.39. P. 58-77 305 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Scalable parallel algorithm for simulation of 3D dynamics of gravitating systems using particle-in-cell method 1 N.V. Snytnikov 1 Institute of Computational Mathematics and Mathematical Geophysics 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). Keywords: gravitating systems, Poisson equation, Vlasov equation, particle-in-cell method, dynamic load balancing. References 1. Hockney, R.W., Eastwood, J.W. Computer Simulation Using Particles. McGraw-Hill. New York. 1981. 540 p. 2. Berezin Ju.A., Vshivkov V.A. Metod chastic v dinamike razrezhennoj plazmy. [Berezin Yu.A. Vshivkov V.A. Particle-in-cell method in the dynamics of low-density plasma] Nauka. Novosibirsk, 1980. 96 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) 306 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 12. Snytnikov N.V., Vshivkov V.A. Metod dekompozitsii oblasti dlya superkomp’yuternogo modelirovaniya gravitiruyushchikh sistem. [Domain decomposition method for supercomuter simulation of gravitating systems] // Russian Supercomputing Days. Proceedings of International Conference. 2015. P. 572-580. 13. Snytnikov V.N. et al. Three-Dimensional Numerical Simulation of a Nonstationary Gravitating N-Body System with Gas // Astronomy Letters. 2004. Vol. 30. P. 124-137. 14. Vshivkov V.A., Snytnikov V.N., Snytnikov N.V. Modelirovanie trekhmernoy dinamiki veshchestva v gravitatsionnom pole na mnogoprotsessornykh EVM [Simulation of three-dimensional dynamics of matter in gravitational field with the use of multiprocessor computer]. // Vychislitelnye tekhnologii (Computational Technologies). 2006. Vol. 11, No. 2. P. 15-27. 15. Samarskii A.A., Andreev V.B. Raznostnye metody dlya ellipticheskih uravneniy [Difference methods for elliptic equations]. Moscow: Nauka, 1976. 352 p. 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, P.24-38. 17. Frigo M., and Johnson S.G. FFTW software. URL: http://www.fftw.org (access date: 01.02.2016) 18. Ayala O., Wang L.P. Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition // Parallel Computing. 2013. Vol.39. P. 58-77. 307