автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Алгоритмы и программы для решения задач динамики многих тел на суперЭВМ

На правах рукописи

КУКШЕВА Эльвира Александровна

АЛГОРИТМЫ И ПРОГРАММЫ ДЛЯ РЕШЕНИЯ ЗАДАЧ ДИНАМИКИ МНОГИХ ТЕЛ НА СУПЕРЭВМ

05.13.18 „Математическое моделирование, численные методы и комплексы программ"

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата технических наук

Новосибирск - 2008

Работа выполнена в Институте катализа им. Г.К. Борескова Сибирского отделения Российской академии наук

Научный руководитель кандидат физико-математических наук,

доцент Сиытников Валерий Николаевич

Официальные оппоненты: доктор физико-математических наук.

профессор Федорук Михаил Петрович

доктор физико-математических наук, доцент Свешников Виктор Митрофанович

Ведущая организация Новосибирский государственный

университет

Защита диссертации состоится " в " А'ОА&ЬЛ- 2008 г. в ~ часов на заседании диссертационного совета Д 003.005.01 при Институте автоматики и электрометрии СО РАН по адресу: 630090, г. Новосибирск, просп. Академика Коптюга, 1.

С диссертацией можно ознакомиться в библиотеке ИАиЭ СО РАН.

Автореферат разослан 3 £>/СГЛЯр&. 2008 г.

Ученый секретарь диссертационного совета, д.ф.-м.н.

Насыров К.А.

Общая характеристика работы

Актуальность темы диссертации.

Проблема расчета траекторий движения многих тел, взаимодействующих между собой, часто встречается в различных областях физики и техники. Эта классическая задача возникает при изучении химических процессов, протекающих в природе, например, в задачах астрокатализа [9]. В этих задачах объектом изучения являются доилаиетные этапы и условия химических превращений вещества, ставшего основой для формирования солнечной системы и органических соединений на поверхности Земли. Одним из главных инструментов этого изучения является математическое моделирование с решением в качестве первого шага задачи динамики многих тел. Математическая модель этой задачи состоит из уравнений Власова-Лиувилля для функции распределения тел по пространству и скоростям и уравнения Пуассона для действующей на тела гравитационной силы. Эта задача решалась методом частица-сетка РМ. В этом методе на введенной сетке по распределению в пространстве большого числа частиц вычисляются значения плотности вещества. Затем, но этим значениям из решения уравнения Пуассона находятся сеточные значения потенциала и соответствующие силы, действующие на тела-частицы. Следом из уравнений движения для частиц определяются их новые координаты. Тем самым, распределение частиц в пространстве изменяется под действием самосогласованного поля, а также из-за других сил. Практически метод частиц представляет собой прямое моделирование динамики многих тел под воздействием собственного и внешнего полей.

Реализация такой модели в трехмерном случае называется ЗБЗУ программой, поскольку используется три пространственные «Б» и три скоростные «V» независимые переменные. Для эксплуатации этих программ с достаточно большим числом частиц, требуются значительные вычислительные ресурсы, поскольку необходимо отслеживать миллионы и миллиарды отдельных траекторий частиц. Размер вычислительной сетки в трехмерных расчетах должен быть не менее 1003 узлов. Поэтому такие программы должны использовать ресурсы суперЭВМ. Но для этого необходимо создать высокоэффективные параллельные алгоритмы и программы, спроектированные так, чтобы в будущем их можно было относительно легко расширять,

При моделировании динамики многих тел РМ-методом возникает проблема постановки граничного условия для уравнения Пуассона. В исходных задачах оно ставится как нулевое значение гравитационного

потенциала на бесконечности. При моделировании в области конечного размера возникает проблема определения значения потенциала на границах при появлении тел-частиц вблизи границ.

Цели работы.

Разработать комплекс программ для решения задач динамики многих тел с применением суперЭВМ. В частности, необходимо было решить следующие задачи:

- разработать метод, позволяющий вычислять потенциал, используя фундаментальное решение уравнения Пуассона;

- реализовать эффективную параллельную процедуру данного метода как в декартовой, так и в полярной системе координат;

- разработать параллельную программу для решения ЗБЗУ задачи многих тел в декартовой системе координат, использующую реализованную процедуру для задания граничных условий.

- использовать эту параллельную процедуру в задачах ЗБ2У для вычисления потенциала на плоскости. Сравнить результаты, полученные этой модифицированной реализацией ЗБ2У задачи с уже имевшимися программами ЗБ2У;

- разработать программу визуализации, ориентированную на представление динамики многих тел.

Научная новизна

работы состоит в следующем:

1. Разработан метод, реализующий дискретный аналог фундаментального решения уравнения Пуассона. Показано, что метод может быть востребован для решения уравнения Пуассона при использовании суперЭВМ с большим, свыше 20, числом процессоров.

2. Разработаны эффективные параллельные реализации этого метода в декартовой и полярной системах координат, которые используют свойство регулярности вычислительной сетки, сокращающее количество операций. Эти реализации использованы в программах ЗБЗУ и ЗБ2У для решения трехмерных задач динамики многих тел.

3. Разработана параллельная программа для решения нестационарных трехмерных задач многих тел ЗБЗУ.

4. Разработанная программа 3D3V использована для практических расчетов и а суперкомпьютерах с делыо поиска квазистациоиарыого решения для диска с вращением. Установлено существование такого численного решения для самогравитирующего бесстолкиовительиого диска и определены его параметры.

Защищаемые положения.

1. Метод решения уравнения Пуассона, который основан на дискретном аналоге его фундаментального решения в декартовой и полярной системах координат, является эффективным при параллельной реализации.

2. Параллельная программа 3D3V для решения нестационарных задач многих тел в трехмерном пространстве.

3. Существует квазистационарное решение задачи о вращении частиц в виде диска, которое получено численно с помощью программы 3D3V на параллельных системах.

Достоверность результатов

подтверждена решением ряда модельных задач, а также сравнением с результатами, полученными другими авторами.

Научная и практическая ценность работы.

Разработанные программы использовались для установления пространственно трехмерных распределений параметров в самогравитирующих эллиптических системах и дисках. Было получено квазистационарное решение для задачи 3D3V. Работа велась в рамках программы Президиума РАН «Происхождение и эволюция звезд и галактик», «Происхождение и эволюция биосферы» (2004-2008 гг.), интеграционного проекта СО РАН №148 (2003-2005 гг.).

Апробация работы.

Основные научные результаты диссертации докладывались и обсуждались на семинарах «Математическое моделирование больших задач» в ИВМиМГ, на Всероссийской конференции «Научный сервис в сети Интернет: решение больших задач» (Новороссийск, 22-27 сентября, 2008), на Всероссийской конференции по вычислительной математике КВМ-2007 (Новосибирск, 18-20 июня, 2007), на 9-th International Conference Parallel Computing Technologies

РаСТ-2007 (Переславль-Залесский, 3-7 сентября, 2007), на международном симпозиуме «Collective Phenomena in Macroscopic Systems» (Como, Italy 4-6 December, 2006) на международном рабочем совещании «Происхождение и эволюция биосферы» (ВОЕ, Новосибирск, 26-29 июня, 2005).

Публикации.

Основное содержание диссертации отражено в 8 работах, список которых приведен в конце автореферата.

Объем и структура диссертации.

Диссертация состоит из введения, шести глав, заключения и списка литературы. Диссертация изложена на 128 страницах, включает библиографический список из 91 наименований работ. Рисунки, формулы и таблицы нумеруются по главам. Библиографические ссылки имеют сквозную нумерацию по всей работе.

Во введении обосновывается актуальность темы диссертации, изложены основные результаты проведенных ранее исследований в области гравитационной физики, сформулирована цель диссертации, приведено краткое содержание работы и полученные результаты.

В первой главе приведен обзор задач гравитационной физики и астрокатализа, а также численные подходы к их решению, которые широко используются в настоящее время.

Вторая глава посвящена исследованию дискретного аналога фундаментального решения уравнения Пуассона с точки зрения его распараллеливания, а также точности по отношению к аналитическим решениям.

Гравитирующие массы распределены в ограниченной плоской области с плотностью /э(т'). Они определяют потенциал, для которого запишем в безразмерных переменных уравнение Пуассона:

Содержание диссертации

Atp — 4-7Г р

Оно имеет фундаментальное решение в виде:

(1)

S

где ¡р(т) - это потенциал, который необходимо вычислить в точке т, |Д| = |т - т'\ - расстояние между точками т и г'. Для дискретного 2В случая в уравнении (2) заменим интеграл суммой и получим:

М1, ч у ^О'Ч

Пусть область решения определяется радиусом Я^ в полярных координатах. Тогда в расчетной области введем равномерную 2Э сетку потенциала [Ф] с шагами кг и 1г„:

Г{ = 1гг ■ г, г =0,1,, Л,, к,. = ^

Ф,,; = Ф(г,;,^) - значение сеточной функции потенциала в узлах сетки

Ф?,,-' находится через суммирование вкладов точечных масс Мц, заданных в центрах ячеек сетки [Ф] но формуле (3). Массы составляют свою сетку [М], сдвинутую относительно сетки потенциала как показано на рисунке 1. Таким

-К

Рис. 1:

образом, Д;,;,;',/ - расстояние от узла (г',/) сетки [Ф] до узла (г,.?) сетки [М]. Поставленная задача определения потенциала в узлах сетки на плоскости прямо связана с методом крупных частиц решения уравнения Власова-Лиувилля [5]. Для задач, включающих решение последнего уравнения, традиционно ориентируются на точность выполнения законов сохранения и расчета полей лучше 1%. Далее будем называть методом (3) численный метод нахождения потенциала с помощью фундаментального решения уравнения Пуассона по формуле (3) на введенных сетках [Ф] и [М].

Для исследования сходимости метода (3) использовались два теста сравнения с аналитическими решениями для бесконечно тонкого кольца и равномерного диска.

Тест с бесконечно тонким кольцом. Аналитическое решение для потенциала кольца определяется формулой:

7Г

ТО I , 1 (4)

и М- = 100, ЛГГ = 200, Мг = 400, Nr = 800,

Л^ = 100 ЛГ„ = 200 ЛГ„ = 400 = 800

0,02 0,5378% 0,1991% 0.1036% 0,0523%

0,5 0.7000% 0.3785% 0,1895% 0,0948%

0,08 1,9130% 0,9992% 0,5022% 0,2501%

1 1,0691% 0,5414% 0,2715% 0,1359%

1,98 0,0001% 0,0002% (1,0001%) 0,0000%

В расчетах кольцо имело радиус Я = 1,01 и массу М = 1. Размерность сетки менялась от (Яг = 100, Л^ = 100), (А/,. = 100, ЛГ, = 200) ,.., до (Ыг = 100,ЛГ^ = 1600). На рисунке 2 показана относительная погрешность метода в процентах для четырех сеток. Как следует из приведенных данных,

Радиус

Рис. 2: Зависимость относительной погрешности от размера сетки для теста с кольцом: (1) - Л'г = 100, Л^ = 100: (2) - к = 100,/V, = 200; (3) - = 100, Л/, = 400. Приведена область наибольших значений относительной погрешности.

есть сходимость к точному решению, и уже на третьей из приведенных сеток точность метода значительно превышает 1% во всей области.

Тест с диском. Аналитическое решение для потенциала равномерного диска определяется формулой:

Ил ?

Г -4- Я-М Г 1 _ , ,

ф г = / -7—Г5\ ■ / I '5

В расчетах радиус диска В,а = 1, его масса М = 1. В таблице представлены значения относительной погрешности, выраженной в процентах, полученного численного решения.

Из данных следует сходимость численного решения к аналитическому решению примерно с первым порядком. Наибольшее отклонение от точного

значения оказалось около края диска. Эта погрешность для сетки ЛГГ = 200, Л^ = 200 меньше 1%, для сетки = 800, Л^ = 800 составляет 0,25%.

Алгоритм распараллеливания. В программной реализации метода (3) задаются массивы потенциала [Ф] и значений точечных масс [М\. Для распараллеливания разрежем массив масс [М] на пр колец и/или секторов, где пр - количество процессоров. В памяти каждого процессора к окажется «его» массив масс М^пр, Лу и полный массив потенциала Ф[ЛГг, Лу. В результате, каждый процессор считает вклад в потенциал тех масс, которые он содержит в памяти по формуле (3). После чего массивы потенциалов из всех процессоров собираются в одном и там суммируются.

Параллельная программа для метода (3) была реализована на языке С II с библиотекой МР1. На рисунке (Рис.3) представлена зависимость ускорения Бр от числа процессоров, где Бр = Т\ - время работы программы на 1-м процессоре, Тр - время работы программы на р процессорах. Ускорение вычислялось на тесте с диском. Из графика видно,

Рис. 3: Зависимость ускорения от числа процессоров для размерностей вычислительных сеток: (1) - Ыг = 100= 100; (2) - Иг = 200, ^ = 200; (Я) - Иг = 500, Л^ = 500

что параллельный алгоритм метода (3) имеет хорошее ускорение на данных тестах для всех представленных размерностей, и при увеличении размерности вычислительной сетки ускорение увеличивается.

В третьей главе диссертации представлена параллельная программа моделирования трехмерной динамики вещества в гравитационном поле, ЗБЗУ.

Система уравнений гравитационной физики без столкновений состоит из кинетического уравнения Власова-Лиувилля и уравнения Пуассона для потенциала гравитационного ноля. Уравнение Власова-Лиувилля в приближении усредненного самосогласованного поля записывается в виде;

ц^+ш.о. №

от ох а V

где /(¿, - зависящая от времени одночастичная функция

распределения по координатам и скоростям, 1? = = -УФ - ускорение частиц единичной массы. Гравитационный потенциал Ф, в котором происходит движение, делится на две части:

Ф = Ф1 + Ф2,

где Ф1 в зависимости от моделируемых условий может представлять потенциал внешних сил, таких как центральное тело, галактическое гало, балдж. Вторая часть потенциала Ф2 определяется совокупным распределением движущихся частиц и удовлетворяет уравнению Пуассона ((? - гравитационная постоянная):

ДФ2 = 4тг Ср, (7)

Для изолированной системы граничное условие состоит в нулевом значении потенциала на бесконечности: Фг->ос = 0.

В некоторый момент времени, который можно принять за начальный, частицы располагаются в расчетной области. Расположение может быть произвольным и, в частности, в виде плоского диска с осесиммстричным распределением поверхностной плотности. Задание начальных скоростей частиц определяет движение диска вокруг центра области. В качестве параметров обезразмеривания выбираются начальный размер гравитирующей системы До, ее масса Мд, и гравитационная постоянная (3.

Отсюда: _

т/ /СМо Во 2

Уо = у иг* = Фо = Ко'

где Ц), ¿о и Фц - характерные величины скорости частиц, времени и потенциала. Все дальнейшие выкладки приведены в безразмерном виде.

Численные методы. Для решения кинетического уравнения Власова используется метод частиц-в-ячейках. В пространственной области в виде параллелепипеда вводится сетка, которая делит область иа ячейки. Модельные частицы, в общем случае с собственной приписанной им массой и другими характеристиками, имеют индивидуальные координаты и могут перемещаться из ячейки в ячейку в соответствии со своими скоростями, определяя функцию распределения по скоростям и координатам. Функция плотности восстанавливается с помощью интерполяции массы каждой частицы в узлы введенной сетки по известным координатам частицы с ядром РЮ. Уравнения движения отдельных частиц имеют вид

й~щ _ ^ <1x1 _ ->

Численное решение этих уравнений осуществляется но следующей схеме, где п - шаг но времени:

У;/'+1 = V," + тРх", Уу"+1 - и,;1 + тЕип, угш = V/1- + г^/',

хп+1 = хп + тух'> + \ у« +1 ^ уп + гщ"+\ ¿» + 1 = + Тг,г»+1

Для нахождения сил, действующих на каждую частицу, используется интерполяция значения сеточной вектор-функции .Р в местоположение частицы с координатами центра ж, у, г, а значения компонент сил, в свою очередь, вычисляются из значений сеточной функции потенциала.

Уравнение Пуассона в прямоугольной сеточной области с регулярной сеткой узлов {х< = Ьхг, у^ = 1гук, г; = г = 0,...,/, к = О,...,Я, I = О,..., Ь], о которой уже упоминалось выше, аппроксимируется схемой второго порядка:

Фи ш - + Ф)-ш , - Щи + +1Л-и , Л.А'.п I - ,1.1+■ 01,1,1-1 , -Л?-+-Л?-+-А?--

Полученная система уравнений решается троекратным применением преобразования Фурье, которое реализуется на основе процедуры быстрого преобразования Фурье. Для этого необходимо независимо найти значение потенциала на границе области.

В соответствии с граничным условием для изолированной системы, состоящем в убывании потенциала от точечной массы как г"1 до нуля на бесконечности, необходимо вычислять значение потенциала тел на границе. Это можно сделать с помощью фундаментального решения уравнения Пуассона. Его дискретный аналог записывается в виде:

= И -Зег (8)

Ш П1',к',11

Значение сеточной функции потенциала находится через суммирование вкладов точечных масс тщы в центрах ячеек введенной сетки. Д-^;/ -расстояние от точечной пцы массы до узла сетки, где вычисляется потенциал

В пашем ЗБ случае нужно вычислить потенциал на шести гранях области моделирования. Поэтому формула (8) дополняется следующими выражениями:

{¿ = 0,..,/, /с = 0,.., К, 1 = 0,..,Ц г' = 0,..,/, к' = 0,..,К, 1' = 0,1' = Ц г' = 0,..,/, к' ~ 0, к' = К, Г = 0, .., Ц г' = 0, г' — к' = 0,..,К, /'= 0,..,1;

ты = Урш> V = ЬхЬуК

- плотность частиц в ячейке, V- объем ячейки. Таким образом, чтобы вычислить потенциал (8) на границе, необходимо определить сеточные

значения плотности и найти сеточные значения потенциала граничных ячеек.

В качестве начальных условий для координат и скоростей частиц могут быть выбраны произвольные распределения. Как пример, рассмотрим одну из наиболее интересных и физически содержательных фигур - модель плоского диска с осесимметричиым распределением поверхностной плотности в виде:

До - начальный радиус диска. В декартовых координатах г = \/{х - Хо)2 + (у - уо)2, г = 2п, где ха, уо, га - центр области. Начальные скорости модельных частиц в диске заданы в соответствии с условием их вращения вокруг центра области. В добавок к этому, в начальный момент времени всем трем компонентам скорости задастся разброс по нормальному закону с дисперсиями с{ьп йм^, (Ь)г.

Параллельная реализация задачи (6),(7) использует 7 трехмерных массивов для хранения сеточных значений функций и 6 одномерных массивов для хранения координат и скоростей порядка 10ь-109 частиц. Параллельный алгоритм состоит в последовательном выполнении следующих основных блоков, каждый из которых в свою очередь выполняется на параллельных процессорах:

1. Инициализация, задание начального распределения частиц и скоростей, вычисление начальной плотности;

2. Вычисление граничных условий для уравнения Пуассона;

3. Решение уравнения Пуассона;

4. Решение уравнения Власова;

5. Вычисление плотности;

6. Сохранение промежуточных результатов.

Начиная со второго, блоки выполняются последовательно в цикле по

времени.

Уравнение Власова. Копии сеточных значений функций потенциала, плотности и трех компонент сил находятся во всех процессорах. Это упрощает распараллеливание, так как позволяет распределять массивы с частицами произвольным образом между процессорами. Недостатком алгоритма является то, что максимально доступный размер расчетных сеток ограничивается объемом оперативной памяти процессорного элемента, например до 256если процессорный элемент имеет 2Гб оперативной памяти.

(9)

При инициализации каждый процессор распределяет в вычислительной области свою часть частиц и вычисляет их начальные скорости. На следующем шаге по времени происходит пересчет положения частиц и их скоростей без обменов массивами частиц между процессорами. Чтобы вычислить плотность, необходимо использовать все частицы, поэтому каждый процессор находит вклад в плотность своих частиц, затем происходит обмен вкладами и их суммирование.

Решение уравнения Пуассона реализовано с помощью библиотеки FFTW, которая позволяет параллельно выполнять многомерные преобразования Фурье для комплексных чисел на многопроцессорных машинах с распределенной памятью и MPI. При распараллеливании выполняется следующее: массив значений потенциала делится но одному из направлений па блоки в соответствии с числом процессоров, каждый процессор выполняет преобразование Фурье для своего блока данных. Затем происходит обмен блоками между процессорами с целью образовать результат во всех процессорах. В итоге, все процессоры имеют массив потенциала с одинаковыми значениями.

Граничные условия. Вычисление потенциала но формуле (8) происходит на 6 граничных плоскостях 3D сетки, но для этого используется весь 3D массив. Оценим трудоемкость вычислений. Для одной плоскости число операций ~ 8п5, где пл- число узлов сетки. Сюда входят операции нахождения расстояния, деления и сложения. Для снижения трудоемкости данной реализации был применен набор приемов. Первый - распараллеливание. Массив плотности разбивается на блоки по числу процессоров N, и каждый процессор вычисляет вклад своей части массива плотности в значение потенциала. Затем происходит пересылка и суммирование вкладов. Пересылаются только те плоскости, для которых выполняются вычисления. Это сокращает число операций примерно в N раз. Используя свойство расчетной сетки (Рис.4(а)), когда ОА = 0\А\ и ОВ = OiBi, можно

о__о,_____

■I "Â

"1 в Ж,

а

Рис. 4:

сэкономить время на вычислении расстояний между узлами сетки. Все эти расстояния можно вычислять при инициализации программы и хранить

\

а

в дополнительном массиве, откуда в процессе счета можно доставать необходимое значение. Этот же массив можно использовать на всех шагах по времени. Число операций при этом, учитывая распараллеливание, сокращается примерно до Для еще большего ускорения на введенную вычислительную сетку можно наложить более грубую сетку, на которой и вычислять потенциал по формуле (8). При этом для вычисления вкладов соседних масс использовать мелкую сетку, а для дальних масс - грубую, как показано на рисунке 4(6), а при вычислении значений па более подробной сетке использовать линейную интерполяцию. В этом случае вместо одного большого массива расстояний размером п'1 хранится два намного меньших массива /с1' и (|)3 , где кл - размерность маленького массива расстояний до ближних масс, а - размерность массива расстояний для дальних масс. Получающееся в результате всех модификаций число операций, можно примерно оцепить как: + ^г. Значение к должно находиться в разумном интервале, с учетом сохранения точности вычислений. В реализованной программе к < у/п, а введенные массивы вместе обычно получаются примерно в 2000 раз меньше основных ЗБ массивов.

Ускорение параллельной программы. Рисунок 5 показывает ускорение параллельной программы на сетке 128'1 и 108 частицами.

зо

Рис. 5: Ускорение параллельной программы ЗОЗУ на сетке 128'' и 108 частицами. Остальные параметры расчета следующие: начальный радиус диска равен 1, его масса - 1, масса центрального тела - 0, дисперсии <йу = 0.5, йЦ, = 0.25, й>г = 0.08, рачмер расчетной области 203. Как следует из рисунка, при числе процессоров более 100 рост-ускорения составляет не более 15% для данных численных параметров.

Четвертая глава посвящена параллельной программе моделирования динамики тел на плоскости. Рассматривается численная модель ЗБ2У задачи, состоящая из уравнения Пуассона для гравитационного поля, уравнения Власова для крупных частиц и газодинамических уравнений. Эта задача уже была решена [5]. Цслыо работы было заменить метод вычисления гравитационного потенциала в этой задаче и сравнить результаты,

полученные обеими программами. В дайной реализации гравитационный потенциал на плоскости находился методом, описанным на странице 6, а в предыдущей реализации он вычислялся комбинированным методом [5]. В результате, были проведены расчеты на суперЭВМ по обеим программам с широким изменением счетных и физических параметров. На основании выполненных исследований можно заключить, что смена вычислительного метода для уравнения Пуассона в задаче 3D2V влияет на детали получаемых результатов. Но в целом поведение системы сохраняется: при одинаковых входных параметрах по обеим программам получаются одинаковые с физической точки зрения результаты.

В пятой главе описывается программа визуализации Gala, специально разработанная для представления результатов моделирования задач гравитационной физики, главным образом в виде фильмов. После небольшого введения, описываются общие принципы работы программы, интерфейс Gala, перечислены доступные графические режимы, описаны допустимые в Gala форматы входных данных. Gala - законченный программный продукт. В настоящий момент она интенсивно используется и автором и коллегами, работающими в данной области.

В шестой главе диссертации представлены результаты численных экспериментов. Для 3D3V задачи приведено квазистационарное решение, эллипсоиды вращения, модельная задача распада гравитирующей среды на сгустки и т.д. Полученное квазистационарное решение системы (6),(7) представлено на рисунке 6. Время Т = 8, начальный радиус диска равен 1, его масса - 1, масса центрального тела - 0, дисперсии dvT = 0.5, dvv = 0.25, dvz = 0.08. Расчет был выполнен на 32 процессорах, размер расчетной области 20'1, размерность вычислительной сетки 128л, число модельных частиц 108.

В заключении перечислены основные результаты, полученные в работе и выводы.

Основные результаты работы.

1. Показано, что реализация метода дискретного аналога фундаментального решения уравнения Пуассона может быть применена для решеиия уравнения Пуассона на параллельных системах с большим, свыше 20, числом процессоров.

2. Разработаны эффективные параллельные реализации этого метода в декартовой и полярной системах координат, которые используют свойство регулярности вычислительной сетки, чтобы сократить

а

б

Рис. ß: Распределение логарифма поверхностной плотности вещества для момента времени Т = 8 (а) - в экваториальной, (б) - в меридиапальной плоскости диска.

количество операций. Эти реализации исиользовапы в программах 3D3V и 3D2V для решения трехмерных задач многих тел.

3. Разработана параллельная программа для решения трехмерных задач многих тел 3D3V, легко переносимая между различными видами кластеров.

4. Проведен ряд расчетов на суперкомпьютерах с целью поиска квазистационарного решения для диска с вращением. В результате установлено существование численного квазистационарного решения для самогравитирующего бесстолкновительного диска.

Основное содержание диссертации опубликовано в работах:

[1] Э.А. Кукшева, В.Н. Снытников Параллельный алгоритм и программа Key7D для решения нестационарных трехмерных задач гравитационной физики. // Вычислительные технологии, Том 12, №1, 2007, с. 35-44

[2] S.Kireev, E.Kuksheva, A.Snytnikov, N.Snytnikov, V.Vshivkov Strategies for Development of a Parallel Program for Protoplanctary Disc Simulation // PaCT 2007, LNCS 4671, pp. 128-139, 2007.

[3] Кукшева Э.А., Снытников В.Н. Параллельная реализация фундаментального решения уравнения Пуассона // Вычислительные технологии, Т. 10, №4, 2005, с. 63-71.

1С

[4] E.A. Kuksheva, V.E. Malyshkin, S.A. Nikitin, A.V. Snytnikov, V.N. Snytnikov, V.A. Vshivkov Supercomputer simulation of self-gravitating media // Future Generation Computer Systems 21 (2005) 749-757

[5] Снытников B.H., Вшивков В.А., Кукшева Э.А., Никитин С.А., неупокоев Е.В., Снытников А.В. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астрономический журнал. 2004, Т.ЗО, №2, с.146-160.

|6] Вшивков В.А., Кукшева Э.А., Никитин С.А., Снытников А.В., Снытников В.Н. О параллельной реализации численной модели физики гравитирующих систем // Автометрия, т. 39, №3, 2003, с.115-123.

[7] Е.A.Kuksheva, V.E.Malyshkin, S.A.Nikitin, A.V.Snytnikov, V.N.Snytnikov, V. A. Vshivkov Numerical Simulation of Self-Organization in Gravitationally Unstable Media on Supercomputers // Lectures in notes. Proc. 7th Int. Conf. PaCT 2003, Nizhni Novgorod, Russia, September 15-19. Springer-Verlag, Berlin, pp. 354-368, 2003

[8] V.N. snytnikov, E.A. Kuksheva Studying Instability of 3D Collision-less Systems on Stochastic Trajectories // COLLECTIVE PHENOMENA IN MACROSCOPIC SYSTEMS Proceedings of the Workshop Villa Olmo, Como, Italy 4 - 6 December 2006.

Цитируемая литература

[9] В.Н. Снытников Абиогенный допланетный синтез пребиотического вещества. // Вестник Российской Академии Наук, Том 77, №3, с. 218-226, 2007,

КУКШЕВА Эльвира Александровна

АЛГОРИТМЫ И ПРОГРАММЫ ДЛЯ РЕШЕНИЯ ЗАДАЧ ДИНАМИКИ МНОГИХ ТЕЛ НА СУПЕРЭВМ

Автореф. дисс. иа соискание ученой степени кандидата технических наук Подписано в печать 01.10.2008. Заказ № 80. Формат 60x84/16. Усл. печ. л. 1,25. Тираж 100 экз. Типография Института катализа им. Г.К. Борескова СО РАН

2007517286

Введение

ГЛАВА 1. Краткий обзор работ по гравитационной физике и астрокатализу.

1.1. Развитие численных подходов.

1.2. Метод частиц-в-ячейках.

1.3. Применение суперЭВМ.

1.4. Работы в области гравитационной физики.

1.5. Численное моделирование задач гравитационной физики.

1.6. Уравнение Пуассона.

1.7. Астрокатализ.

ГЛАВА 2. Метод определения потенциала, основанный на фундаментальном решении уравнения Пуассона.

2.1. Исходные уравнения.

2.2. Сходимость и точность метода.

2.2.1. Тест с бесконечно тонким кольцом.

2.2.2. Тест с диском.

2.2.3. Сравнение с другими методами.

2.3. Параллельная реализация.

2.3.1. Алгоритм распараллеливания.

2.3.2. Эффективность распараллеливания.

2.4. Обсуждение результатов.

ГЛАВА 3. Параллельная реализация численной модели гравитирующей изолированной системы, задача 3D3V.

3.1. Постановка задачи.

3.2. Численные методы.

3.2.1. Уравнение Власова.

3.2.2. Уравнение Пуассона.

3.2.3. Граничные условия для уравнения Пуассона.

3.2.4. Начальные условия.

3.3. Параллельная реализация.

3.3.1. Уравнение Власова.

3.3.2. Решение уравнения Пуассона.

3.3.3. Граничные условия.

3.3.4. Ускорение параллельной программы.

3.4. Тестирование параллельной реализации.

3.4.1. Сохранение энергии.

3.5. Программная реализация.

3.6. Выводы.

ГЛАВА 4. Использование фундаментального решения для задачи динамики тел на плоскости. Квазитрехмерный случай.

4.1. Исходная система уравнений.

4.2. Численные методы.

4.2.1. Уравнение Власова-Лиувилля.

4.2.2. Газодинамические уравнения.

4.2.3. Уравнение Пуассона.

4.3. Параллельная реализация.

4.3.1. Распараллеливание уравнение Пуассона.

4.3.2. Распараллеливание метода частиц.

4.3.3. Газовая динамика.

4.3.4. Схема работы программы.

4.4. Оценка эффективности параллельной реализации.

4.5. Влияние счетных параметров.

4.6. Сравнение двух реализаций задачи 3D2V.

4.6.1. «Кольца».

4.6.2. «Спирали».

4.7. Выводы.

ГЛАВА 5. Программа визуализации.

5.1. Как работает программа Gala.

5.2. Интерфейс Gala.

5.2.1. Выбор файлов.

5.2.2. Файл конфигурации.

5.2.3. Создание анимации. Запуск визуализации.

5.2.4. Помощь и подсказки.

5.3. Графические режимы.

5.4. Файл конфигурации.

5.4.1. Карта.

5.4.2. Вектора.

5.4.3. Поверхность.

5.5. Форматы файлов данных.

5.5.1. Текстовый формат задания функции зависимости от 2-х переменных.

5.5.2. Бинарный формат задания функции зависимости от 2-х переменных.

5.5.3. Бинарный формат задания функции зависимости от 3-х переменных.

Проблема расчета траекторий движения многих тел, взаимодействующих между собой, часто встречается в различных областях физики и техники. Эта классическая задача возникает при изучении химических процессов, протекающих в природе, например, в задачах астрокатализа [72, 61]. В этих задачах объектом изучения являются допланетные этапы и условия химических превращений вещества, ставшего основой для формирования солнечной системы н органических соединений на поверхности Земли. Одним из главных инструментов этого изучения является математическое моделирование с решением в качестве первого шага задачи динамики многих тел. Математическая модель этой задачи состоит из уравнений Власова-Лиувилля для функции распределения тел по пространству и скоростям и уравнения Пуассона для действующей на тела гравитационной силы. Эта задача решалась методом частица-сетка РМ [5, 6]. В этом методе на введенной сетке по распределению в пространстве большого числа частиц вычисляются значения плотности вещества. Затем, по этим значениям из решения уравнения Пуассона находятся сеточные значения потенциала и соответствующие силы, действующие на тела-частицы [6]. Следом из уравнений движения для частиц определяются их новые координаты. Тем самым, распределение частиц в пространстве изменяется под действием самосогласованного поля, а также из-за других сил. Практически метод частиц представляет собой прямое моделирование динамики многих тел под воздействием собственного и внешнего полей.

Реализация такой модели в трехмерном случае называется 3D3V программой, поскольку используется три пространственные «D» и три скоростные «V» независимые переменные. Для эксплуатации этих программ с достаточно большим числом частиц, требуются значительные вычислительные ресурсы, поскольку необходимо отслеживать миллионы и миллиарды отдельных траекторий частиц. Размер вычислительной сетки в трехмерных расчетах должен быть не менее 1003 узлов. Поэтому такие программы должны использовать ресурсы суперЭВМ. Но для этого необходимо создать высокоэффективные параллельные алгоритмы pi программы, спроектированные так, чтобы в будущем их можно было относительно легко расширять.

При моделировании динамики многих тел РМ-методом возникает проблема постановки граничного условия для уравнения Пуассона. В исходных задачах оно ставится как нулевое значение гравитационного потенциала на бесконечности. При моделировании в области конечного размера возникает проблема определения значения потенциала на границах при появлении гравитирующих тел вблизи границ.

Цель работы заключается в разработке комплекса программ для решения задач динамики многих тел с применением суперЭВМ. В частности, необходимо было решить следующие задачи:

- разработать метод, позволяющий вычислять потенциал, используя фундаментальное решение уравнения Пуассона;

- реализовать эффективную параллельную процедуру данного метода как в декартовой, так и в полярной системе координат;

- разработать параллельную программу для решения 3D3V задачи многих тел в декартовой системе координат, использующую реализованную процедуру для задания граничных условий.

- использовать эту параллельную процедуру в задачах 3D2V для вычисления потенциала на плоскости. Сравнить результаты, полученные этой модифицированной реализацией 3D2V задачи с уже имевшимися программами 3D2V;

- разработать программу визуализации, ориентированную на представление динамики многих тел.

Научная новизна работы состоит в следующем:

1. Разработан метод, реализующий дискретный аналог фундаментального решения уравнения Пуассона. Показано, что метод может быть востребован для решения уравнения Пуассона при использовании суперЭВМ с большим, свыше 20, числом процессоров.

2. Разработаны эффективные параллельные реализации этого метода в декартовой и полярной системах координат, которые используют свойство регулярности вычислительной сетки, сокращающее количество операций. Эти реализации использованы в программах 3D3V и 3D2V для решения трехмерных задач динамики многих тел.

3. Разработана параллельная программа для решения нестационарных трехмерных задач многих тел 3D3V.

4. Разработанная программа 3D3V использована для практических расчетов на суперкомпьютерах с целью поиска квазистационарного решения для диска с вращением. Установлено существование такого численного решения для самогравитирующего бесстолкновительного диска и определены его параметры.

Защищаемые положения:

1. Метод решения уравнения Пуассона, который основан на дискретном аналоге его фундаментального решения в декартовой и полярной системах координат, является эффективным при параллельной реализации.

2. Параллельная программа 3D3V для решения нестационарных задач многих тел в трехмерном пространстве.

3. Существует квазистационарное решение задачи о вращении частиц в виде диска, которое получено численно с помощью программы 3D3V на параллельных системах.

Достоверность результатов подтверждена решением ряда модельных задач, а также сравнением с результатами, полученными другими авторами.

Научная и практическая ценность работы.

Разработанные программы использовались для установления пространственно трехмерных распределений параметров в самогравитирующих эллиптических системах и дисках. Было получено квазистационарное решение для задачи 3D3V. Работа велась в рамках программы. Президиума РАН «Происхождение и эволюция звезд и галактик», «Происхождение и эволюция биосферы» (2004-2008 гг.), интеграционного проекта СО РАН №148 (2003-2005 гг.).

Структура и объем работы.

Диссертация состоит из введения, шести глав, заключения и списка литературы. Диссертация изложена на 126 страницах, включает библиографический список из 91 наименований работ. Рисунки, формулы и таблицы нумеруются по главам. Библиографические ссылки имеют сквозную нумерацию по всей работе.

4.7. Выводы.

Создана параллельная реализация задачи 3D2V, в которой уравнение Пуассона вычисляется методом, основанным на его фундаментальном решении.

Проведены расчеты на суперЭВМ по обеим программам с широким изменением счетных и физических параметров.

В результате выполненных исследований можно заключить, что смена вычислительного метода для уравнения Пуассона в задаче 3D2V влияет на детали получаемых результатов. Но в целом поведение системы сохраняется: при одинаковых входных параметрах по обеим программам получаются одинаковые с физической точки зрения результаты.

BICSSRAS V5m*nik<.v. VVihK^ov. EKukihM о BIC 5B PAS V Str^n&ov. WshMtov. E Kukshevd

T .1 В BIC SB RA5 V.5rrXr*oV. V.VjhiUtov, £ Kukrfmu о BIC 50 RAS V Snvtmlfoi:, V.vjhvltov. £ K(*iheve

Рис. 4.6. Расчет «спирали». Распределение плотности пыли в зависимости от времени (по горизонтали) и программы (по вертикали).

ГЛАВА 5 ПРОГРАММА ВИЗУАЛИЗАЦИИ.

Работа над моделированием динамики гравитирующих систем ведется уже достаточно давно [2, 59, 13]. Нужно отметить, что проведение численных экспериментов на суперЭВМ обладает некоторой спецификой, особенно, если речь идет о нестационарных задачах. Прежде чем получить заветную «картинку» приходится выполнить массу различных манипуляций: удаленно запустить приложение, удаленно получить файлы-результаты (часто это сотни), которые нужно заархивировать и «скачать» с удаленного кластера. Затем нужно настроить программу визуализации для своих данных. Иногда требуется переформатировать данные и т.д. Все это необходимо проделывать много раз. Поэтому, в какой-то момент возникла потребность максимально автоматизировать этот процесс, создав программное обеспечение. Желательно было визуализировать нестационарные процессы не в виде отдельных статических изображений, а в виде фильмов. Был проведен анализ доступных на тот момент систем визуализации и не было найдено той, которая в точности подходила бы под наши нужды, а более менее подходящие программы были достаточно дорогими. Решено было создать собственную программу визуализации для представления результатов моделирования задач гравитационной физики. Она должна быть простой в использовании, уметь быстро работать с большими объемами информации (сотни и тысячи Мб), уметь производить фильмы в распространенных форматах. Это и было реализовано в программе Gala. Нужно подчеркнуть, что Gala - это межплатформенная реализация, то есть ее можно откомпилировать и под платформу Windows и под X Windows Unix-систем, потому что Gala реализована на межплатформенной библиотеке Qt.

Gala интенсивно используется в настоящий момент и автором и коллегами, работающими в данной области. Изображения полученные Gala использовались в оформлении данной работы, в опубликованных статьях ([13], [14]), в презентациях конференций, отчетах, на семинарах.

5.1. Как работает программа Gala.

Суть работы программы: на вход подаются файлы данных, а результатом являются соответствующие изображения в графическом формате. Входные данные могут быть как в текстовом формате, так и в бинарном. Выходные изображения могут быть в различных графических форматах. Форма выходного изображения также может быть различной и определяется понятием «графический режим». Для задания всех этих параметров работы программы существует файл конфигурации. Анимация в Gala создается в формате MNG и MPEG2 (MNG - это более новый, свободный в распространении и более экономичный чем GIF формат, нацеленный на использование в интернете).

5.2. Интерфейс Gala.

После запуска программы пользователь видит типичное окно приложения с меню и окно редактирования входного списка файлов.

Все важные функции дублированы в панели инструментов. Пункт меню Windows/File list и Windows/Config editor позволяет переключаться между окном редактирования файла конфигурации и окном редактирования входного списка файлов. Ж Shew

Name j Ext | Size | Date | Time | Path |

File Edit Windows Help

Рис. 5.1

5.2.1. Выбор файлов. Если входные файлы в нужном формате, то можно начать работу с выбора файлов, пункт меню File/Open directory или Ctrl-О или значок панели инструментов. Выберите в появившемся окне директорию с Вашими файлами, и все файлы из этой директории окажутся в списке файлов. Далее Вы можете удалять ненужные файлы из списка с помощью Edit/Del selected files или Ctrl-X, выделяя файлы мышью или SHIFT-стрелки. Отменить выделение можно с помощью Edit/Unselect all или Ctrl-U. В полученный список можно добавить другие директории аналогичным образом через Edit/Add, directories или Ctrl-А.

Можно открыть директорию рекурсивно с помощью меню File/Open directory recursively, в результате в списке окажутся все поддиректории и файлы поддиректорий, выбранной директории.

Список файлов можно отсортировать по: имени, расширению, размеру, дате, времени выбрав соответствующий заголовок мышью. Если в список попадут «неправильные файлы», то при обработке эти файлы будут пропущены.

5.2.2. Файл конфигурации. После выбора директории файлов данных, необходимо выбрать файл конфигурации. Предположим, мы уже имеем отредактированный файл, тогда его просто открываем через File/Load config или Ctrl-L или значок

После этого он появляется в нашем редакторе, и если какие-то настройки в нем пользователя не устраивают, его можно отредактировать и сохранить изменения либо File/Save config либо Ctrl-S. Отредактированный файл можно сохранить с новым именем File/Save config as. Если нужно создать новый файл конфигурации, то выбираем File/New config file/N, где N графический режим». Будет создана «рыба» файла, которую затем можно будет изменить под нужды пользователя.

Очередность операций выбора файлов данных и файла конфигурации не имеет значения.

5.2.3. Создание анимации. Запуск визуализации. Когда файлы данных в списке уже имеются, файл конфигурации тоже в наличии и отредактирован, можно приступать к визуализации (созданию анимации). Но перед этим пользователь может задать имя будущему файлу анимации, если конечно, он намерен получить именно анимацию, а не просто набор графических файлов. Имя задается через File/Save animation as. Если имя не задать, то оно будет равно имени первого файла данных + ".mng". Запуск анимации/визуализации осуществляется через File/Create animation или Ctrl-R или кнопкой панели инструментов. Если Вы хотите видеть как протекает процесс, то необходимо нажать кнопку

Это позволит видеть окно с формирующимися изображениями. Если Вы не желаете видеть процесс формирования изображений, то кнопка должна быть отжатой. После работы, в соответствующих директориях появятся графические файлы в соответствии с выбранными опциями заданного файла конфигурации.

5.2.4. Помощь и подсказки. По мере возможности все кнопки программы снабжены выпадающими подсказками желтого цвета. Работает кнопка «что это». Пункты меню Help/about F1 вызывает окно с документацией пользователя, Help/What's this равносилен кнопке «что это».

5.3. Графические режимы.

Графическим режимом (ГР) здесь называется то, в каком виде мы видим наши данные. Gala имеет б ГР, то есть умеет шестью способами отображать входные данные. Некоторые из них представлены на рисунке 5.2

5.4. Файл конфигурации.

Каждый ГР имеет параметры, изменяющие изображение, которое будет генерироваться. Эти параметры задаются в файле конфигурации. Файл конфигурации - это текстовый файл с древовидной, XML-подобной формой записи параметров. Он задает параметры для работы визуализатора программы Gala. Без указания файла параметров программа ничего не делает. Для разных графических режимов необходимы разные файлы конфигурации. Обычно файл конфигурации состоит из нескольких логических разделов: общие, связанные с входными данными, с выходными данными и с конкретным графическим режимом. Получить «рыбу» файла конфигурации можно в самой Gala через пункт меню File/New config file. Рассмотрим файлы конфигурации для некоторых графических режимов.

5.4.1. Карта. В таблице 5.4.1. приведены примеры файлов конфигурации для карт, векторов и поверхностей. Пример файла конфигурации для получения изображения, показанного на рисунке 5.2(а,) находится в первой колонке таблицы.

Поясним значения этих строк. <global> - начало тэга с общей информации; д

Рис. 5.2. Графические режимы программы Gala.

ЗАКЛЮЧЕНИЕ

В работе исследован дискретный аналог фундаментального решения уравнения Пуассона с точки зрения эффективности его распараллеливания, а также точности по отношению к аналитическим решениям. Показано, что метод может быть востребован для решения уравнения Пуассона при использовании суперЭВМ с числом процессоров свыше 20, особенно для расчета граничных условий в полностью трехмерных задачах. Метод идеально подходит для работы именно на параллельных системах.

Разработана параллельная программа 3D3V для решения трехмерных задач гравитационной физики без столкновений для изолированных систем, которая показала свою эффективность на существующих суперЭВМ. При вычислении гравитационного потенциала, и в этом особенность данной реализации, граничные условия находились с помощью дискретного аналога фундаментального решения уравнения Пуассона.

Для процедуры вычисления гравитационного потенциала на границе был разработан экономичный по времени параллельный алгоритм, использующий регулярность вычислительной сетки, чтобы сократить количество операций. Его применение при вычислении потенциала на границе дает возможность устанавливать границы вычислительной области ближе к моделируемой системе. При этом не нарушается условие убывания потенциала до нуля на бесконечности.

Разработана параллельная программа для решения задач гравитационной физики 3D2V. В этой реализации гравитационный потенциал на плоскости находится с помощью дискретного аналога фундаментального решения уравнения Пуассона, что позволяет уменьшать расчетную область и не находить граничные условия. Эта программа использовалась для сравнения и верификации с результатами других реализаций [2, 3, 13]. Выяснилось, что результаты полученные всеми программами не противоречат друг другу.

Была создана программа визуализации Gala для представления результатов моделирования задач гравитационной физики. Она специализирована на представлении именно таких данных, умеет быстро работать с большими объемами информации (сотни и тысячи Мб), умеет производить фильмы в распространенных форматах. Gala - это межплатформенная реализация, то есть ее можно откомпилировать и под платформу Windows и под X Windows Unix-систем.

Проведен ряд расчетов на параллельных системах, в результате чего было получено численное, квазистационарное решение для самогравитирующего бесстолкновительного диска с вращением. Таким образом, созданный комплекс программ может служить основой для дальнейшей разработки параллельных программ численного моделирования задач динамики многих тел.

Автор выражает благодарность В.А. Вшивкову за внимание и поддержку.

1. Снытников В.Н., Вшивков В.А., Кукшева Э.А., Никитин С.А., неупокоев Е.В., Снытников А.В. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астрономический журнал. 2004, Т.ЗО, №2, с.146-160.

2. В.Н. Снытников, В.Н. Пармон, В.А. Вшивков, Г. И. Дудникова, С.А. Никитин, А.В. Снытников. Численное моделирование гравитационных систем многих тел с газом / / Вычислительные технологии. 2002, Т. 7, №3, с. 72-84.

3. Вшивков В.А., Кукшева Э.А., Никитин С.А., Снытников А.В., Снытников В.Н. О параллельной реализации численной модели физики гравитирующих систем // Автометрия, 2003, Т.39, №3, с.115-123.

4. Б. Олдер, С. Фернвах и М. Ротенберг Вычислительные методы в физике плазмы. М.:Мир, 1974. с.186-187.

5. Хокни Рю, иствуд Дж. Численное моделирование методом частиц. М.:Мир, 1987. 639 с.

6. Григорьев Ю.Н., Вшивков В.А. Численные методы "частицы-в-ячейках". Новосибирск: Наука. Сибирская издательская фирма РАН, 2000. 184 с.7. поляченко В.Л., фридман Л.М. Равновесие и устойчивость гравитирующих систем. Мир.: Наука, 1976. 446 с.

7. Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике.М: Наука, 1982, 392 с.

8. Вшивков В.А., Снытников В.Н., Снытников Н.В. Моделирование трехмерной динамики вещества в гравитационном поле на многопроцессорных ЭВМ. // Вычислительные технологии, 2006, Т. 11, №2, с. 15-27.

9. Кукшева Э.А., Снытников В.Н. Параллельная реализация фундаментального решения уравнения Пуассона // Вычислительные технологии, 2005, Т. 10, Ж, с. 63-71.

10. BORIS J.P. Relativistic plasma simulation optimization of a hybrid code. // Proc. of the Fourt Conf. Num. Sim. Plasmas. Washington, 1970. P.3.

11. 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, 21, 749-757

12. Кукшева Э.А., Снытников В.Н. Параллельный алгоритм и программа Key7D для решения нестационарных трехмерных задач гравитационной физики. // Вычислительные технологии, 2007, Том 12, №1, с. 35-44

13. Официальный сайт FFTW http://www.fftw.org/

14. Aarseth S., Lecar М. Computer simulations of stellar systems. Ann. Rev. Astron. Astrophys., 1975, 13, 1.

15. РОЖАНСКИЙ И.Д. Анаксагор. М.: Наука, 1972.

16. ДЖИНС Дж. Движение миров. М.: Гостехиздат, 1933.

17. SPITZER L. The Dynamical Evolution of Globular Clusters. // Princeton: Princeton Univ. Press., 1987, (Имеется перевод: Спитцер JI. Динамическая эволюция шаровых скоплений. М.: Мир, 1990)

18. Miller R.H., Prendergast К.Н. Stellar Dynamics in a Discrete Phase Space. // Astrophys J. 1968, vol. 151, pp. 699-709.24. hohl F., R.W. hockney A Computer Model of Disks of Stars. // J. Comput. Phys. 1969, vol. 4, pp. 306-323.

19. Miller R.H. Validity of disk galaxy simulations. // J. Сотр. Phys., 1976, 21, 400-437.

20. SELLWOOD, J.A. Bar instability and rotation curves. // Astron. Astrophys., 1981, 99, 362-374.

21. BROWNRIGG D.R.K. Computer Modelling of Spiral Structure in Galaxies. PhD thesis, 1975, Univ. Reading, Engl.

22. SELLWOOD J.A. The Art of N-body Building. // Ann. Rev. Astron. Astrophys. 1987, 25, 151.32. страуструп Б. Язык программирования с++. М.: Невский диалект, 1999, с. 991

23. Harlow F.H., Dickman D.O., Harris D.E., Martin R.E. Two diman-sional hydrodynamic calculations. // Los Alamos Scie. Lab. Rep. NLA-2301., 1959.34. харлоу ф.х. Численный метод «частиц-в-ячейках» для задач гидродинамики. М.: Мир, 1967, 383 с.

24. Теоретические основы и конструирование численных алгоритмов задач математической физики. Под ред. К.И. Бабенко М.: Наука. 1979, 295 с.

25. ЧАН Теренс Системное программирование на С++ для UNIX. Пер. с англ.-К.: Наука. Издательская группа BHV, 1999, 592 с.41. miller R.H., prendergast К.Н. Stellar Dynamics in a Discrete Phase Space. // Astrophys J., 1968, vol. 151, pp. 699-709.

26. Miller R.H., Prendergast K.H. and W.J. Quirk Numerical Experiments on Spiral Structure. // Astrophys. J., 1970, vol. 161, pp. 903-916.

27. Flanders P.A4., D.J. Hunt, S.F. Reddaway and D. Parkinson Efficient High Speed Computing with the Distributed Array Processor, in «High Speed Computer and Algorithm Design». // Academic Press Ltd., London,1977, pp. 113-128.

28. Miller R.H. Free Collapse of a Rotating Sphere of Stars. // Astrophys J.,1978, vol. 223, pp. 122-128

29. HOCKNEY R.W. Super-Computer Architecture, in «Infotech State of the Art Report: Future Systems I» Infotech International, Maidenhead, Berks, 1977, pp. 277-305.

30. Поляченко В.JI., Фридман A.M. Равновесие и устойчивость гравитирующих систем. Мир.: Наука, 1976, 446 с.

31. ФРИДМАН A.M. Предсказание и открытие новых структур в спиральных галактиках. // Успехи физических наук, 2007, Т. 177, №2, с. 121-148.48. bertin G. Dynamics of Galaxies. // Cambridge University Press, 2000.

32. MacFarland Т., Couchman H.M.P., Pearce F.R., Pilchmeier J. A. New Parallel P3M Code for Very Large-Scale Cosmological Simulations. // New Astron. 1998, V. 3, №8, pp. 687-705

33. Bode P., Ostriker J., Xu G. The Tree-Particle-Mesh N-body Gravity Solver. // Astrophys.J. Suppl. 2000, 128, 561-570

34. Lanzoni В., Guiderdoni В., Mamon G.A., Devriendt J., Hatton J. GALICS VI Modelling Hierarchical Galaxy Formation in Clusters. // Mon. Not. R. Astron. Soc., 2005, 1-17.

35. Sellwood J. Galaxy dynamics by N-body simulation. // Rutgers University preprint, 1996.

36. Bureau M., athanassoula E. Bar Diagnostics in Edge-On Spiral Galaxies. III. N-body Simulations of Disks. // Astrophys J., 2005, vol. 626, pp. 159-173.

37. Demmel J.W. Applied Numerical Linear Algebra. SIAM. 1997.

38. Снытников B.H., Пармон В.Н. Жизнь создает планеты? // Наука из первых рук. 2003, Т.0, с.2-15.

39. Kireev S., Kuksheva е., Snytnikov A., Snytnikov N., Vshivkov V. Strategies for Development of a Parallel Program for Protoplanetary Disc Simulation // PaCT 2007, 2007. LNCS 4671, pp. 128-139

40. Snytnikov V.N., Kuksheva E.A. Studying Instability of 3D Collision-less Systems on Stochastic Trajectories. // COLLECTIVE PHENOMENA IN MACROSCOPIC SYSTEMS. Proceedings of the Workshop. 4-6 Dec. 2006, Villa Olmo Como, Italy

41. Снытников А.В., Снытников H.B., Кукшева Э.А. Математическое моделирование на суперЭВМ ранних этапов эволюции планет. // International Workshop Biosphere Origin and Evolution June 26-29, 2005, Novosibirsk, Russia

42. Snytnikov V.N., Vshivkov V.A., Parmon V.N. Solar Nebula as a Global Reactor for synthesis of Prebiotic Molecules // 11th Int. Conf. on the Origin of Life. July 5-12, 1996. p.65, Orleans, France,

43. Вшивков В.А., Никитин С.А., Снытников В.Н. Исследование неустойчивостей бесстолкновительных систем по стохастическим траекториям // Письма в ЖЭТФ. 2003, Т.38, №6

44. Вланшет Ж., Саммерфилд М. Qt3: программирование GUI на С++. КУДИЦ-Образ, 2005, 464 с.

45. Сайт библиотеки OpenGL (русскоязычный) http://opengl.org.ru/

46. Sellwood J.A., Merrit David Instabilities of counterrotating stellar disks. // The Astrophysical Journal. 1994, 430-550,

47. Sellwood J.A., merritt d. Bending Instabilities of Stellar Systems. // The Astrophysical Journal. 1994, 551-567,

48. Shigeru Wakita, Minoru Sekiya Numerical Simulations of the Gravitational Instability in the Dust Layer of a Protoplanetary Disk using a Thin Disk Model. // Astrophysical Journal. 2008, Vol. 675, Issue 2, pp. 1559-1575

49. XAKEH Г. Информация и самоорганизация. Макроскопический подход к сложным системам. Перевод с английского. Серия "Синергетика: от прошлого к будущему". Изд.2, доп. М.: УРСС, 2005, 248 с.

50. Шикин А. В., боресков А.В. Компьютерная графика. Динамика, реалистические изображения. М.: ДИАЛОГ-МИФИ, 1995, 288 с.