CC BY
42
Вычислительные технологии
Том 12, № 1, 2007
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ И ПРОГРАММА Key7D ДЛЯ РЕШЕНИЯ НЕСТАЦИОНАРНЫХ ТРЕХМЕРНЫХ ЗАДАЧ ГРАВИТАЦИОННОЙ ФИЗИКИ*
Э.А. Кукшева, В.Н. Снытников Институт катализа им. Г. К. Борескова СО РАН, Новосибирск, Россия e-mail: kuksheva@catalysis.nsk.su, snyt@catalysis.nsk.su
An effective parallel Key7D code for solution of 3D problems of collisionless gravitational physics in Cartesian coordinates is developed. The discrete analogue of fundamental solution of Poisson's Equation is used to determine boundary conditions during calculation of the gravitational potential. For this procedure an efficient parallel algorithm is developed. A quasistationary solution of the system is presented using the Key7D program.
Введение
Численное моделирование широко применяется для решения задач бесстолкновительной гравитационной физики. К ним можно отнести, в частности, исследования коллективной динамики тел в околозвездных дисках и нестационарных процессов в галактиках, изучение галактических скоплений и крупномасштабной структуры Вселенной. Основную роль в таком моделировании играют методы частиц [1-3]. Среди них важное место принадлежит методу частица-сетка PM. В этом методе на введенной сетке по распределению в пространстве большого числа частиц вычисляется плотность вещества. Затем из решения уравнения Пуассона находятся сеточные значения потенциала и соответствующие силы, действующие на частицы [3]. Из уравнений движения для частиц определяются их новые координаты. Тем самым распределение частиц в пространстве изменяется под действием самосогласованного поля, а также других сил. Практически метод частиц предполагает математическое моделирование из физически первых принципов.
Вычислительные эксперименты с использованием программ по PM-методу позволяют находить детальные немаксвелловские функции распределения частиц в шестимерном фазовом пространстве, отслеживать тонкие эффекты с участием резонансных частиц и определять с контролируемой точностью физические характеристики моделируемых бес-столкновительных систем. Одно из важнейших приложений метода связано с изучением динамики развития гравитационной неустойчивости в различных физических системах.
* Работа выполнена при финансовой поддержке программы Президиума РАН "Происхождение и эволюция звезд и галактик", Министерства образования и науки РФ (гранты № РНП.2.1.1.1969 и № НШ 6526.2006.3).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
Из дисперсионного анализа волн в неограниченных средах с гравитацией следуют устойчивость гравитационных возмущений с длинами волн меньше джинсовской длины и их неустойчивость при больших масштабах, если такие возмущения могут развиваться [4]. В расчетах физических систем с проявлениями гравитационной неустойчивости размер сеточного шага должен находиться в области устойчивости. Поэтому он должен быть значительно меньше характерной длины физической системы, в данном случае — джинсовской длины Ад. Если предположить, что воспроизводимые длины волн находятся в интервале 0.1Ад ... 10Ад, то размер вычислительной сетки в трехмерных 3Б по пространственным переменным расчетах должен быть не менее 1003 узлов. Для обеспечения уровня флукту-аций плотности в 1 % в ячейке должно быть до 104 частиц. Тогда общее количество частиц для указанной сетки оценивается в 1010. В принципе современные суперкомпьютеры имеют технические ресурсы для решения задач с указанными параметрами и воспроизведением детальных функций распределений в многомерных фазовых пространствах. Но для этого необходимо создать высокоэффективные параллельные алгоритмы и программы.
В работе [5] описана параллельная программа для решения нестационарных трехмерных по пространственным переменным задач в цилиндрических координатах. Эта программа предназначена для проведения вычислительных экспериментов, в которых важно сохранение азимутальной симметрии или ее нарушение, к примеру, в галактических или околозвездных дисках. Но ряд других задач, в частности о взаимодействиях галактик или двойных звезд с дисками, предпочтительнее решать в декартовых координатах. Однако для этих задач возникает проблема постановки граничных условий.
В исходной задаче для изолированной системы граничные условия ставятся как нулевое значение потенциала на бесконечности. При использовании области конечного размера с существованием гравитирующих тел вблизи границ нужно дополнительно определить на границах значение потенциала. В работе [6] нами предложено использовать фундаментальное решение уравнения Пуассона для нахождения потенциала многих тел на плоскости в так называемом 3Б2У случае, когда движение вдоль одного из направлений из рассмотрения исключается. Там же показано, что разработанная параллельная реализация перспективна для проведения вычислительных экспериментов на существующих суперкомпьютерах, причем число модельных частиц может достигать порядка 1010. Используя эту идею в случае декартовой сетки для определения граничных условий и необходимые компоненты программного обеспечения современных суперкомпьютеров, возможно создать полностью трехмерную численную модель нестационарной гравитационной физики без столкновений, которую условно будем в дальнейшем называть 3Б3У задачей. Таким образом, целью настоящей работы является создание и изучение параллельного алгоритма 3Б3У задачи гравитационной физики в декартовой системе координат с заданием потенциала на границе для изолированной системы.
1. Постановка задачи
Система уравнений гравитационной физики без учета столкновений состоит из кинетического уравнения Власова — Лиувилля и уравнения Пуассона для потенциала гравитационного поля [7]. Уравнение Власова — Лиувилля в приближении усредненного самосогласованного поля записывается в виде
д£ д£ д£ _ 0
дЬ дх ОН '
где f — зависящая от времени одночастичная функция распределения по координа-
там и скоростям; а = Р/ш = —VФ — ускорение частиц единичной массы. Гравитационный потенциал Ф, в котором происходит движение, делится на две части:
Ф = Ф1 + Ф2,
где Ф1 в зависимости от моделируемых условий может представлять потенциал внешних сил, таких как центральное тело, галактическое гало, балдж. Вторая часть потенциала Ф2 определяется совокупным распределением движущихся частиц и удовлетворяет уравнению Пуассона (О — гравитационная постоянная)
ДФ2 = 4пОр,
которое в декартовой системе координат запишется в виде
дФ. + + = (2)
дх2 ду2 дг2 '
Для изолированной системы граничное условие состоит в нулевом значении потенциала на бесконечности = 0).
В некоторый момент времени, который можно принять за начальный, частицы располагаются в расчетной области. Расположение может быть произвольным, в частности в виде плоского диска с осесимметричным распределением поверхностной плоскости. Задание начальных скоростей частиц определяет движение диска вокруг центра области. В качестве параметров обезразмеривания выбираются размер гравитирующей системы Д°, ее масса М° и гравитационная постоянная О. Отсюда
V, = ,/ ОМ0, *0 = ^' Ф° = у02-
Здесь V,, и Ф° — характерные величины скорости частиц, времени и потенциала. Все дальнейшие выкладки приведены в безразмерном виде.
2. Численные методы
Уравнение Власова. Для решения кинетического уравнения Власова используется метод частиц-в-ячейках [3]. В пространственной области в виде параллелепипеда вводится сетка, которая делит область на ячейки. Модельные частицы, в общем случае с собственной приписанной им массой и другими характеристиками, имеют индивидуальные координаты и могут перемещаться из ячейки в ячейку в соответствии со своими скоростями, определяя функцию распределения по скоростям и координатам. Функция плотности восстанавливается с помощью интерполяции массы каждой частицы в узлы введенной сетки по известным координатам частицы с ядром PIC. Уравнения движения отдельных частиц имеют вид
difi F dcci dt mi dt i
Численное решение этих уравнений осуществляется по следующей схеме, где n — шаг по времени:
vxn+1 = vxn + rFxn, vyn+1 = vyn + rFyn, vzn+l = vzn + rFzn,
хп+1 = хп + ТУХП+1, уп+1 = уп + т^п+1, гп+1 = гп + т^п+1.
Для нахождения сил, действующих на каждую частицу, используется интерполяция значения сеточной вектор-функции Г в местоположение частицы с координатами центра х, у, г, а значения компонент сил, в свою очередь, вычисляются из значений сеточной функции потенциала.
Уравнение Пуассона. Уравнение Пуассона в прямоугольной сеточной области с регулярной сеткой узлов (хг = Нхг, ук = 'ук, гг = кгI, г = 0,...,1, к = 0,...,К, I = 0,...,Ь}, о которой уже упоминалось выше, аппроксимируется схемой второго порядка [3]:
фг+1,к,1 — 2фг,к,1 + фг-1,к,1 . фг,к+1,1 — 2фг,к,1 + фг,к-1,1 . фг,к,1+1 — 2фг,к,1 + фг,к,1-1 л
-Г~2- + -Г~2- + -Г~2- = 4ПРг,к,1.
Их "у 'г
Полученная система уравнений решается применением троекратного преобразования Фурье, которое реализуется на основе процедуры быстрого преобразования Фурье. Для этого необходимо независимо найти значение потенциала на границе области.
Граничные условия для уравнения Пуассона. В соответствии с граничным условием для изолированной системы, состоящим в убывании потенциала от точечной массы как г-1 до нуля на бесконечности, необходимо вычислять значение потенциала тел на границе. Это можно сделать с помощью фундаментального решения уравнения Пуассона. Его дискретный аналог записывается в виде
фг',к',1' =^ „Ш . (3)
г,к,г пг',к',1'
Значение сеточной функции потенциала фг'кк',г' находится через суммирование вкладов
г,г,кЛ
точечных масс тг,к,г в центрах ячеек введенной сетки. Щ к' у — расстояние от точечной тг,к,г массы до узла сетки, где вычисляется потенциал фг',к',г'. В нашем 3Б случае нужно вычислить потенциал на шести гранях области моделирования. Поэтому формула (3) дополняется следующими выражениями:
г = 0,..,1, к = 0,..,К, I = 0,..,Ь,
г = о,..,1, к' = о,..,к, V = 0, V = ь,
г' = 0,..,1, к' = 0, к' = К, V = о,..,ь, г = 0, г = I, к' = 0,..,к, V = 0,..,ь,
тг,к,г = Ург,к,г, У = 'х'у кх,
где Ргкк,,\ — плотность частиц в ячейке; У — объем ячейки. Таким образом, чтобы вычислить потенциал на границе (3), необходимо определить сеточные значения плотности и найти сеточные значения потенциала граничных ячеек. Как показано в [6], этот способ вычисления потенциала на границе позволяет проводить расчеты при нахождении тел вблизи и на границах расчетной области.
Начальные условия. В качестве начальных условий для координат и скоростей частиц могут быть выбраны произвольные распределения. В качестве примера рассмотрим одну из наиболее интересных и физически содержательных фигур — модель плоского диска с осесимметричным распределением поверхностной плотности в виде
г
-«иг ' Г; £ (4)
Здесь Д0 — начальный радиус диска. Значения цилиндрических координат частиц по (р и г задаются в соответствии с
р € тив [0; 2п], г = г0,
где тив — функция, генерирующая случайные числа на интервале [0; 2п]; г0 — центральная плоскость области по г. Начальные скорости модельных частиц в диске заданы по нормальному закону со значением дисперсии по всем трем направлениям пространства — вьг, вь^, вьх.
3. Параллельная реализация
Параллельная реализация задачи (1), (2) использует семь трехмерных массивов для хранения сеточных значений функций и шесть одномерных массивов для хранения координат и скоростей порядка 106 ... 109 частиц. Параллельный алгоритм состоит из следующих основных блоков, каждый из которых выполняется на параллельных процессорах:
1) инициализация, задание начального распределения частиц и скоростей, вычисление начальной плотности;
2) вычисление граничных условий для уравнения Пуассона;
3) решение уравнения Пуассона;
4) решение уравнения Власова;
5) вычисление плотности;
6) сохранение промежуточных результатов.
Начиная со второго, блоки выполняются последовательно в цикле по временной координате.
Уравнение Власова. Представленный параллельный алгоритм не использует декомпозицию области. Копии сеточных значений функций потенциала, плотности и трех компонент сил находятся во всех процессорах. Это упрощает распараллеливание, так как позволяет массивы с частицами распределять произвольным образом между процессорами. Недостатком алгоритма является то, что максимально доступный размер расчетных сеток ограничивается объемом оперативной памяти процессорного элемента, например до 2563, если процессорный элемент имеет 2 Гбайта оперативной памяти. При инициализации каждый процессор распределяет свою часть частиц и вычисляет их начальные скорости. Далее происходит пересчет положения частиц и их скоростей на следующем шаге по времени без обменов массивами частиц между процессорами. Чтобы вычислять плотности, необходимо использовать все частицы, поэтому каждый процессор находит вклад своих частиц в плотность, затем происходят обмен и суммирование этих вкладов.
Решение уравнения Пуассона реализовано с помощью библиотеки FFTW, которая позволяет параллельно выполнять многомерные преобразования Фурье для комплексных чисел на многопроцессорных машинах с распределенной памятью и MPI [8]. Средствами FFTW массив значений потенциала делится по одному из направлений на блоки в соответствии с числом процессоров, каждый процессор выполняет преобразование Фурье для своего блока данных. Затем осуществляется обмен блоками между процессорами с целью сбора результатов во всех процессорах.
Граничные условия. Вычисление потенциала по формуле (3) выполняется на шести граничных плоскостях 3D сетки, но для этого используется весь 3D массив. Оценим
трудоемкость вычислений. Для одной плоскости число операций ~ 8п5, где п3 — число узлов сетки. Вычисления включают операции нахождения расстояния, деления и сложения. Для снижения трудоемкости данной реализации был применен набор приемов. Первый — распараллеливание. Массив плотности разбивается на блоки по числу процессоров N, и каждый процессор вычисляет вклад своей части массива плотности в значение потенциала. Затем происходят пересылка и суммирование вкладов. Пересылаются только те плоскости, для которых выполняются вычисления. Это сокращает число операций примерно в N раз. Используя свойство расчетной сетки (рис. 1, а), когда О А = 0\и О В = 0\В\, можно сэкономить на вычислении расстояний между узлами сетки. Все эти расстояния вычисляются при инициализации программы и хранятся в дополнительном массиве. Затем в процессе счета оно берется из массива расстояний. Этот массив используется на всех шагах по времени. Число операций при этом, учитывая распараллеливание, сокращается примерно до 3п5 /N. Для еще большего ускорения на введенную вычислительную сетку накладывается более грубая сетка. На ней и вычисляется потенциал по формуле (3). Для вычисления вкладов соседних масс используется мелкая сетка, а для дальних масс — грубая, как показано на рис. 1, б, а при вычислении значений на более подробной сетке используется линейная интерполяция. В этом случае вместо одного большого массива расстояний размером п3 хранится два намного меньших массива к3 и (п/к)3 , где к3 — размерность маленького массива расстояний до ближних масс, а (п/к)3 — размерность массива расстояний для дальних масс. Получающееся в результате всех модификаций
Рис. 1.
Рис. 2. Ускорение параллельной программы в зависимости от числа процессоров.
9п5 2п2 0
число операций можно примерно оценить как + ——. Значение к должно находить-
^к2 к2
ся в разумном интервале с учетом сохранения точности вычислений. Обычно к < у/п. Введенные массивы вместе обычно получаются примерно в 2000 раз меньше основных 3Б массивов.
Ускорение параллельной программы. Ускорение параллельной программы на сетке 1283 и 108 частиц показано на рис. 2. Как из него следует, при числе процессоров более ста рост ускорения составляет не более 15 % для данных численных параметров.
4. Тестирование параллельной реализации
Параллельная программа верифицировалась на следующих тестах.
1. Решения уравнения Пуассона и уравнения Власова тестировались отдельно друг от друга на аналитических решениях.
2. В расчетах контролировались импульс, момент импульса и полная энергия, которая сохраняется лучше 1 %.
3. Сходимость решения проверялась на последовательности сгущающихся сеток. На рис. 3 показана поверхностная плотность в центральном сечении по диагонали диска — оси У. На рис. 3, а применена логарифмическая шкала, он иллюстрирует сходимость при измельчении шага сетки при 108 частиц и то, как при измельчении вычислительной сетки становятся явными более тонкие структуры решения по амплитуде в центре. Сходимость решения при увеличении числа частиц приведена на рисунке 3, б. Сдвиг пиков плотности в зависимости от числа частиц происходит в пределах одной сеточной ячейки, что присуще методу частиц в ячейках РМ.
4. Сравнение результатов расчетов с результатами аналогичной программы, описанной в [5]. Обе программы были запущены с одинаковыми начальными параметрами, вплоть до полного совпадения начальных значений координат и скоростей частиц на сетке 643 с 106 частицами. На рис. 4 показана поверхностная плотность в центральном сечении по У для времени Т = 4. Несмотря на то что для решения системы (1), (2) применялись раз-
Рис. 3. Зависимость распределения поверхностной плотности в центральном сечении по У от счетных параметров: а — зависимость логарифма поверхностной плотности от вычислительной сетки, сетка 643, 1283 и 2563 (кривые 1-3 соответственно), число частиц 108; б — зависимость поверхностной плотности от числа частиц, 106, 108 и 4 • 108 (кривые 1-3 соответственно); сетка 1283; время Т = 2.
Рис. 4. Распределение поверхностных плотностей в центральном сечении по диагональной оси У, полученное с применением разных программ: кривая 1 — программы, описанной в данной статье, кривая 2 — программы, описанной в [5]. Параметры расчетов: сетка 643, 106 частиц, время Т = 4.
личные численные методы и системы координат, полученные трехмерные нестационарные решения различаются между собой лишь в пределах расчетной декартовой ячейки.
5. Квазистационарное решение
Созданная программа Key7D использовалась для поиска квазистационарных решений системы (1)-(3). Эти расчеты производились на кластере (BladeCenter) в Межведомственном суперкомпьютерном центре РАН и Сибирском суперкомпьютерном центре. Приведем результат одного из таких расчетов, выполненных на кластере BladeCenter, имеющем порядка 900 процессоров IBM PowerPC-970+ c оперативной памятью 2 Гбайта на процессорный элемент и сеть Myrinet. На рис. 5 представлена поверхностная плотность в экваториальной и меридиональной плоскостях диска на времени T = 8, начальный радиус диска равен 1, его масса 1, масса центрального тела 0, дисперсии dvr = 0.5, dv^ = 0.25, dvz = 0.08. Расчет выполнен на 32 процессорах, размер расчетной области 203, размерность вычислительной сетки 1283, число модельных частиц 108. Расчет выполнялся в течение 15.5 ч.
б
Рис. 5. Распределение логарифма поверхностной плотности вещества для момента времени Т = 8
в экваториальной (а) и меридиональной (б) плоскостях диска.
Рис. 6. Поведение полной энергии системы во времени.
На рис. 6 показано сохранение полной энергии системы. Изменение относительной погрешности составило менее 1 %. При указанных начальных данных система вращается в самосогласованном гравитационном поле с угловой скоростью, зависящей от расстояния до оси вращения. К моменту времени Т _ 12 система выполнила свыше двух оборотов, рассчитанных по угловой скорости среды в точке с поверхностной плотностью 0.1 от аналогичной плотности в центре системы (см. рис. 3, а). Таким образом, получен пример квазистационарного бесстолкновительного самогравитирующего диска без предположений о существовании дополнительных гравитационных потенциалов (темной материи, балджа, гало и т. д.).
Заключение
Описан разработанный параллельный код решения трехмерных задач гравитационной физики без столкновений для изолированных систем, который показал свою эффективность на существующих суперЭВМ. При вычислении гравитационного потенциала, и в этом особенность данной реализации, граничные условия находились с помощью дискретного аналога фундаментального решения уравнения Пуассона. Для этой процедуры разработан экономичный по времени параллельный алгоритм, использующий для сокращения количества операций свойство регулярности вычислительной сетки. Его применение для вычисления потенциала на границе дает возможность устанавливать практически любые размеры области. При этом не нарушается условие убывания потенциала до нуля на бесконечности. Проведен ряд расчетов на суперкомпьютерах с целью поиска квазистационарного решения для диска с вращением. В результате установлено существование численного решения для самогравитирующего бесстолкновительного диска.
Авторы выражают благодарность В.А. Вшивкову за полезные консультации по созданию алгоритма метода частиц.
Список литературы
[1] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987. 639 с.
[2] Miocchi P., Capuzzo-Dolcetta R. An efficient parallel tree-code for the simulation of self-gravitating systems // Astronomy & Astrophysics. 2002. Vol. 382, N 2. P. 758-767.
[3] Григорьев Ю.Н., Вшивков В.А. Численные методы "частицы-в-ячейках". Новосибирск: Наука. Сиб. изд. фирма РАН, 2000. 184 с.
[4] ПолячЕнко В.Л., Фридман Л.М. Равновесие и устойчивость гравитирующих систем. М.: Мир; Наука, 1976. 446 с.
[5] Вшивков В.А., Снытников В.Н., Снытников Н.В. Моделирование трехмерной динамики вещества в гравитационном поле на многопроцессорных ЭВМ // Вычисл. технологии. 2006. Т. 11, № 2. С. 15-27.
[6] Кукшева Э.А., Снытников В.Н. Параллельная реализация фундаментального решения уравнения Пуассона // Вычисл. технологии. 2005. Т. 10, № 4. С. 63-71.
[7] Снытников В.Н., Вшивков В.А., Кукшева Э.А., Никитин С.А. и др. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астроном. журн. 2003. Т. 29, № 12.
[8] Официальный сайт FFTW http://www.fftw.org/
Поступила в редакцию 26 июля 2006 г., в переработанном виде —15 декабря 2006 г.
