Содержание к диссертации
Введение
2 Обзор литературы 9
2.1 Обзор исследований астрофизических дисков 9
2.1.1 Описание структуры галактики 9
2.1.2 Основные понятия теории волн плотности 10
2.1.3 Оценка влияния столкповительных процессов 12
2.1.4 Оценка релятивистской поправки 13
2.1.5 Протопланетные диски 14
2.2 Обзор вычислительных методов 16
2.2.1 Обзор методов и программ моделирования гравитируюших систем 16
2.2.2 Обзор методов решения уравнения Пуассона 18
2.2.3 Обзор работ по многосеточному методу 19
3 Исходные уравнения 21
3.1 Уравнение Власова-Лиувилля 22
3-2 Уравнение Пуассона 22
3-3 Уравнения газовой динамики 23
4 Численные методы 24
4.1 Решение кинетического уравнения 24
4.2 Решение уравнений газовой динамики 25
4.2.1 Уравнения Эйлерова этапа метода крупных частиц 26
4.2.2 Разностная схема для эйлерова этапа метрда крупных частиц 27
4.2.3 Уравнения для Лагранжева этапа метода крупных частиц . 30
4.2.4 Решение уравнения температуры газа 34
4.2.5 Тесты 38
4.3 Решение уравнения Пуассона 43
4.3.1 Граничные условия 43
4<3.2 Переход к системе линейных алгебраических уравнений * 45
4.3.3 Методы решения 49
4.3.4 О преимуществе итерационных методов 50
4.3.5 Метод блочной последовательной верхней релаксации . 51
4.3.6 GMRES 53
4.3.7 Многосеточный метод 54
4.3.8 Тесты 58
5 Параллельная реализация численных методов 59
5.1 Применение сборочной технологии 59
5.2 Декомпозиция области 60
5.3 Метод частиц 62
5.4 Уравнение Пуассона 63
5.5 Оценка времени работы программы 65
6 Численные эксперименты 67
6.1 Устойчивость диска G7
6.1.1 Инкременты неустойчивости 67
6.1.2 Эволюция диска 70
6.1.3 Влияние динамической температуры на устойчивость диска 73
6.1.4 Модифицированная скорость Тумре 75
6.1.5 Разогрев диска 76
6.1.6 Сравнение с наблюдениями 77
6.2 Неустойчивость двухфазного диска 78
6.2Л Описание экспериментов 78
6.2.2 Азимутальная нсустойчипость 79
6.2.3 Порог неустойчивости 81
6.3 СолитоиоподоОные структуры 83
6.3.1 Параметры расчетов 83
6.3.2 Влияние сетки па качество расчета 84
6.3.3 Поведение солитоноподобных волн 87
6.3.4 Структура солитоноподобной волны 90
7 Заключение
- Основные понятия теории волн плотности
- Уравнение Пуассона
- Разностная схема для эйлерова этапа метрда крупных частиц
- Декомпозиция области
Введение к работе
Диски являются фундаментальным понятием астрофизики, причина этого состоит в том, что существенная часть наблюдаемого вещества находится в системах, по внешнему виду напоминающих диски [1], к ним относится большинство галактик. Но и в самих галактиках диски встречаются достаточно часто: наблюдаются протозвезды - совсем молодые, только рождающиеся звезды, окруженные дисками. Из газопылсвых дисков возникают планеты, дискообразные системы имеются вокруг планет - примерами могут служить кольца Сатурна, Юпитера, Нептуна. Итак, в природе имеется целая иерархия дисков. По-видимому, впервые понятие астрофизического диска возникло в предложенной Лапласом теории образования Солнечной системы.
Физические условия в перечисленных выше системах чрезвычайно многообразны, Вещество дисков выступает в форме звезд, газовых облаков, газа, пыли; спектр температур лежит в широчайших пределах от десятков до миллионов градусов. Дискообразные объекты объединяет между собой прежде всего то, что существуют они благодаря совместному действию вращения и гравитации. Диски различаются не только пространственными масштабами (от сотен километров до сотен килопарсек), но и характерными временами процессов -от долей секунд до миллиардов лет. Несмотря па все бросающиеся в глаза различия, диски во многих случаях обнаруживают большое сходство в отношении протекающих в них физических процессов.
В физике дисков существует несколько нерешенных проблем, в частности:
проблема формирования планет из газопылевого протоплаиетиого диска, в особенности газовых гигантов [2]
проблема эволюции вещества в галактиках [3]
проблема происхождения спирального узора в галактических дисках [4]
Для решения этих проблем необходимо численное моделирование [5].
Вышеперечисленные проблемы имеют общую особенность: в дисках галактик, также как и в протоплапетных дисках, можно выделить две фазы вещества. В протопланетном диске это газ и твердые тела размером порядка метра, в галактическом диске - газ и звезды. В динамике газовой фазы существенную роль играют столкновения. Другая фаза, состоящая из звезд либо твердых тел, которая будут далее условно называться фазой частиц, является практически бесстолкновителыюй. Таким образом, актуальным является моделирование двухфазного диска с гравитацией.
Задача динамики диска является нестационарной и существенно трехмег>-ной- Как показано в [1], при наличии центрального тела, масса которого более чем в 5 раз превышает массу диска, диск, содержащий газовую фазу и фазу
частиц, может рассматриваться как двумерный. Для решения этой задачи была создана гибридная модель [б, 7], состоящая из расчета бесстолкновительного движения частиц путем решения уравнения Власова-Лиувилля и решения уравнений газовой динамики в гравитационном поле. Гравитационное поле в этой модели складывается из двух частей: самосогласованное поле диска, которое определяется из уравнения Пуассона и поле центрального тела.
Сложность этой модели заключается в принципиальной пространственной трехмерности уравнения Пуассона и необходимостью проследить за индивидуальным движением большого числа частиц. Количество модельных частиц определяет точность восстановления функции распределения по скоростям и по координатам, содержащейся в уравнении Власова—Лиувилля. Для персональных ЭВМ число отслеживаемых частиц достигает миллионов во всей области с приблизительно 10* ячейками, что обеспечивает в среднем сотни частиц в ячейках и уровень флуктуации, к примеру, плотности до десяти процентов, так как флуктуации зависят от числа частиц как jV~1|/2. Для снижения уровня флуктуации макроскопических параметров хотя бы к 0,1% количество частиц в расчетной области необходимо увеличивать на четыре порядка. Таким образом, для проведения расчетов необходимо использование суперЭВМ.
Цель работы заключается в создании параллельного алгоритма для моделирования двухфазных дисковых гравитирующих систем и численной реализации этого алгоритма. Для этого необходимы:
разработка быстрых параллельных методов решения уравнения Пуассона;
параллельная реализация метода частиц в ячейках;
» модификация метода крупных частиц для корректного учета границы газ-вакуум;
исследование численной модели;
проведение вычислительных экспериментов на суперЭВМ.
Двухфазный диск с гравитационным полем может считаться двумерным лишь в определенном диапазоне физических параметров (прежде всего - массивное центральное тело) [1], Расчеты за пределами этого диапазона проводятся с целью исследования свойств самой численной модели.
Научная новизна работы состоит в следующем:
создание параллельного алгоритма для моделирования двухфазных дисковых гравитирующих систем и его численная реализация;
построение метода решения уравнения Пуассона, оптимизированного для решения нестационарной задачи. Основной особенностью предложенного метода являются возможность ускорения счета за счет использования потенциала, полученного на предыдущем временном шаге;
разработка существенно параллельного метода решения уравнения Пуассона на основе преобразования Фурье;
увеличение ускорения параллельной программы в результате применения многосеточного метода;
обнаружение в вычислительных экспериментах неустойчивости, обусловленной взаимным влиянием двух фаз в диске.
Научная ценность работы заключается п создании комплекса программ для моделирования двухфазных дисковых гравитирующих систем на суперЭВМ, исследовании свойств численной модели и изучении устойчивости двухфазных дисковых гравитирующих систем при различных значениях физических параметров. Расчеты динамики двухфазного диска в двумерной постановке могут быть полезны как необходимый фактический материал при переходе к трехмерной постановке задачи, существенно расширяющей возможности численного моделирования. Полученная реализация численной модели двухфазных дисковых гравитирующих систем на суперЭВМ может быть использована при моделировании протопланетпых или галактических дисков.
Представленные в диссертации исследования проводились в соответствии с планами ИВМиМГ СО РАН по программе суперЭВМ, в рамках программы 25-2 Президиума РАН "Происхождение и эволюция биосферы", при поддержке интеграционного проекта СО РАН помер 148 и федеральной целевой программы "Интеграция"; поддерживались Российским Фондом Фундаментальных Исследований (проект 02-01-00864, 05-01-00665), а также при частичной поддержке БАСА, грант JURRJSS NAG 5-8717.
Достоверность результатов подтверждается тестированием, которое заключается в следующем. По отдельности метод решения уравнения Пуассона и газодинамическая схема были протестированы на стационарных задачах, имеющих аналитическое решение. Сходимость численных методов проверена на последовательности измельчающихся сеток. Дополнительно правильность расчета контролируется с помощью проверки выполнения законов сохранения энергии, массы, момента импульса. Также было произведено тестирование на модельных задачах, близких по физической постановке к изучаемым явлениям и допускающим аналитическое решение. Результаты вычислительных экспериментов по моделированию устойчивости звездной подсистемы сравнивались с аналитическим критерием Тумре [8]. Решение уравнений газовой динамики проверено на задаче о распаде разрыва и истечении газа в вакуум.
Вычислительная часть диссертационной работы отдельно опубликована в центральных реферируемых журналах из списка ВАК и доложена на ряде конференций. Результаты проведенного исследования докладывались на различных (в том числе и международных) конференциях и семинарах отдельно по организации и по проведению расчетов {обоснование расчетной части работы)
и отдельно по анализу полученных физических результатов. Результаты работы достаточно полно опубликованы в рецензируемых центральных научных журналах по списку ВАК соответствующего профиля. Автор защищает:
параллельный алгоритм для моделирования двухфазных дисковых грави-тирующих систем и его численную реализацию;
параллельный метод решения уравнения Пуассона;
метод решения уравнения Пуассона, оптимизированный для решения нестационарной задачи;
модификацию метода крупных частиц решения задач газовой динамики.
Структура и объем работы.
Диссертация состоит из введения, трех глав, заключения и списка литературы. Диссертация изложена па 107 страницах, включает библиографический список из 110 наименований работ. Рисунки, формулы, таблицы и библиографические ссылки имеют сквозную нумерацию но всей работе.
Во введении обоснована актуальность темы диссертации, изложены основные результаты проведенных ранее исследований по моделированию и наблюдениям протопланетных дисков, содержится обзор проблем численного решения задач, рассматриваемых в диссертации, сформулирована цель диссертации и приведено краткое изложение полученных результатов. Дается краткое описание дисковых галактик и протопланетных дисков, приводится краткое изложение основных положений теории воли плотности, которые необходимы для анализа протекающих в диске процессов. Также дано обоснование применимости классической динамики для моделирования диска с гравитацией и бесстолкно-вителыюго приближения для расчета движения пылевых частиц»
В разделе "Исходные уравнения" выписаны исходные уравнение - уравнение Пуассона, уравнения газовой динамики и уравнение Власова-Лиувилля,
В первой главе описывается построение численных методов решения исходных уравнений.
В разделе 4.1 согласно монографии [9] описывается решение уравнения Власова-Лиувилля методом частиц в ячейках. Описан переход от кинетического уравнения к уравнениям движения модельных частиц, приводятся уравнения движения в цилиндрической системе координат с поправкой Бориса [10]. Подчеркивается внутренний параллелизм метода.
В разделе 4.2 описывается реализация метода крупных частиц Белоцерков-ского-Давыдова применительно к динамике газа з плоскости диска. Приведены уравнения эйлерова и лагранжева этапов метода, описала численная схема решения уравнений на эйлеровом этапе. Реализована схема коррекции энергии газа с использованием уравнения температуры. Проведено тестирование работы
разностной схемы на задаче о распаде разрыва, а также на задаче роста неустойчивости с заданной длиной волны.
В разделе 4.3 разрабатывается метод решения уравнения Пуассона для плоского протопланетного диска в цилиндрической геометрии. Предложены два комбинированных метода: быстрое преобразование Фурье в сочетании с методом последовательной верхней релаксации и быстрое преобразование Фурье в сочетании с многосеточным методом. Основной особенностью предложенных комбинированных методов являются возможность учета ускорения решения за счет использования потенциала, полученного на предыдущем временном шаге. Вследствие использования быстрого преобразования Фурье методы обладают внутренним параллелизмом, который дает возможность минимизировать объем межпроцессорных коммуникаций. Представлена методика расчета граничных условий для уравнения Пуассона, оптимальным по соотношению точность - производительность является метод, использующий моменты инерции диска. Указаны преимущества использования метода быстрого преобразования Фурье; внутренний параллелизм и сокращение времени расчета. Реализованы различные методы решения плохо обусловленной системы линейных алгебраических уравнении, возникающей при дискретизации уравнения Пуассона, Приведены результаты сравнения, наиболее быстрым оказывается многосеточный метод.
Во второй главе описывается параллельная реализация разработанных численных алгоритмов. Параллельная программа строится на основе сборочной технологии параллельного программирования (раздел 5.1), Рассмотрены ла-грапжева и эйлерова декомпозиция расчетной области (раздел 5.2), Реализовано два варианта лагранжевой декомпозиции: разбиение па подобласти в физическом пространстве по радиальной координате и в пространстве азимутальных фурьс-гармоник по волновому числу. В силу того, что движение частиц считается двумерным, удалось получить более простую и эффективную реализацию метода частиц (раздел 5.3).
Третья глава посвящена вычислительным экспериментам по моделированию динамики двухфазного диска. Особое внимание уделяется вопросам устойчивости диска (раздел ОД)- Приведен различные критерий устойчивости, разработаны средства для визуализации инкрементов неустойчивости. Проведен ряд экспериментов с однофазным (состоящим только из частиц) и двухфазным дисками. Приведены полученные в ходе экспериментов распределения плотности пыли и газа, инкременты неустойчивости, проведен анализ соответствия теории и данным наблюдений.
В разделе G.2 описаны вычислительные эксперименты по моделированию двухфазного диска, в которых была получена неустойчивость, обусловленная взаимным влиянием двух фаз. Отмечено, что каждая из фаз по отдельности являются устойчивыми. Проанализированы возможные причины ее возникновения, отмечено сходство с двухпотоковой неустойчивостью в плазме.
В разделе 6.3 рассмотрены получештые в вычислительных экспериментах
уединенные волны плотности, похожие по виду на солитоны. Далее эти волны будут называться солитоноподобными структурами или солитоноподобными волнами. Описаны эффекты, связанные с ними: поглощение солитоноподобной волны и полная неподвижность этой волны в потоке вещества. Таким образом показано, что солитоноподобные структуры являются именно волнами ПЛОТНО' сти, а не сгустками вещества. Рассмотрена структура солитоноподобной волны: приведены графики плотности и скорости частиц и распределение плотности газа.
В заключении представлены основные результаты, полученные d работе и выводы.
Основные результаты работы:
Создан параллельный алгоритм для моделирования двухфазных дисковых гравитирующих систем и его численная реализация.
Разработан метод решения уравнения Пуассона, оптимизированный для решения нестационарной задачи;
Построена модификация метода крупных частиц решения задач газовой динамики;
Разработан существенно параллельный метод решения уравнения Пуассона для моделирования динамики двухфазного диска с гравитационным полем;
Исследованы свойства солитоноподобных структур в двухфазном диске с гравитационным полем;
В вычислительных экспериментах обнаружена неустойчивость, обусловленная взаимным влиянием двух фаз в диске с гравитационным полем.
Личный вклад соискателя заключается в обсуждении постановок задач, разработке адекватных численных алгоритмов и методов решения, создании и тестировании алгоритмов и программ, проведении расчетов, интерпретации результатов численного моделирования. Все выносимые на защиту результаты принадлежат лично автору. Представление изложенных в диссертации и выносимых на защиту результатов, полученных в совместных исследованиях, согласовано с соавторами.
Апробация работы- Основные научные результаты диссертации опубликованы в централыюй печати [7], [11] - [18], докладывались и обсуждались на семинаре "Математическое и архитектурное обеспечение параллельных вычислений" под руководством д.т.и. В.Э.Малышкина, на семинаре "Математическое моделирование больших задач" под руководством д.ф.-млі. В,А.Вшивкова, на Сибирской школе-семинаре по параллельным вычислениям (Томск, ноябрь
2003), на школе-семинаре ''Распределенные и кластерные вычисления" (Красноярск, 17-19 сентября 2002), на международном научно-практическом семинаре "Высокопроизводительные вычисления на кластерных системах" (Нижний Новгород, 20-24 ноября 2001,13-15 сентября 2003) па Международной Конференции по Вычислительной Математике (Новосибирск, 24-28 июня 2002, 21-25 июня 2004), на международной конференции Parallel Computing Technologies (Нижний Новгород, 15-19 сентября 2003, Красноярск, 5-9 сентября 2005), а также на международном рабочем совещании "Происхождение и эволюция биосферы" (26-29 июня 2005), Кроме того, были сделаны доклады на конференциях молодых ученых, ИВТ СО РАН (25-26 декабря 2000), ИЯФ СО РАН (апрель 2002):ИМВиМГСОРАН (март 2002, апрель 2005) и на международной научной студенческой конференции "Студент и научно-технический прогресс" (Новосибирск, 2000, 2001, 2002). Материалы докладов опубликованы в сборниках трудов конференций [19] - [23].
Автор выражает глубокую признательность Валерию Николаевичу Снытни-кову за постановку задачи, Виталию Андреевичу Вшивкову за осуществление научного руководства и поддержку па всех этапах работы над диссертацией и Виктору Эммануловичу Малышкииу за ряд цепных замечаний и консультаций.
2 Обзор литературы
Основные понятия теории волн плотности
Впервые идея о волнах плотности как о механизме образования спиральных рукавов в галактиках была высказана Линдбладом [25]. Аналогично волнам пространственного заряда п плазме, малые возмущения плотности рассматриваются в виде
Е Si (г) exp [І{ІОІ — ттыр)} здесь г, р - полярные координаты, m - число спиральных рукавоп, OJ - частота волн, которая может быть представлена в виде ш = mftp — 7 здесь Іїр - угловая скорость спирального узора, 7 инкремент.
Минимальная длина, на которой может развиваться гравитационная неустойчивость - это так называемая джинсовская длина качественно - это минимальный размер капли, который удерживается вместе вследствие самогравитации. В плазме аналогичной величиной является дебаев-СКИЙ радиус TD = 4тгпе
Дисперсионное уравнение для волн плотности проще всего получить в гидродинамическом приближении. (ш - mil)2 = к2 + k2c2 - 1vG r\\k\\ (1) где с - тепловая скорость звезд, к - эпициклическая частота / ldlnu к-2и\ /14-"—— V 2d\nr
Это уравнение дает хорошее качественное описание поведения волн, в частности, оно показывает два механизма стабилизации диска: вследствие высокой динамической температуры или вращения. Однако уравнение (1) применимо только в приближении больших длин волн.
\_( 1)ЯШ)
В работе Лина и Illy [2G] была впервые показана возможность существования спиральных волн, получено дисперсионное соотношение позволяющее по известным параметрам галактики и скорости волны воспроизвести геометрию спирального узора. fcj 1 (2) X Sl("/n)2-lJ здесь 1п - модифицированная функция Бесселя порядка n kj - критическое волновое число, соответвующее джинсовской длине kj — тр, х = г-" и V безразмерная частота it) — 771 V = AJ1 смысл которой в следующем: это частота, с которой вещество галактики попадает в рукава спирального узора.
Также в [26] были указаны особенности воли, вызванные неоднородностью системы (резонаисы Линдблада, коротациониый резонанс). При коротации звезды и волна вращаются с одинаковой угловой скоростью. При внутреннем (внешнем) резонансе Линдблада угловая скорость звезд больше (меньше), чем угловая скорость волны; однако звезды все время проходят пик волны в той же самой фазе и оказывают когерентное воздействие на нее - это явление аналогично затуханию Ландау [27] в плазме. Условие резонансов имеет вид [3]: ilv = ш 4- s— где s — 0 при коротаний, s — 1 и -1 при внешнем и внутреннем резонансе Линдблада соответственно. Видно, что при этом знаменатель в правой части уравнения (2) обращается в 0. и - угловая скорость звезд в диске, 1Р - угловая скорость волны.
В задаче моделирования генерации волн плотности диск галактики может считаться бесстолкновительным. Время, в течение которого столкновениями можно пренебречь, пропорционально числу частиц. В данном разделе это будет показано на примере диска галактики [3], для которой время релаксации сопоставимо с возрастом вселенной.
Время релаксации по энергии TJJ определяется из соотношения & = TD J(Avp)22ir(vnn)pdp то есть это время, за которое поперечная кинетическая энергия частицы mVp/2 становится сравнимой с изначальной кинетической энергией mv /2. Отсюда получается формула, аналогичная случаю рассеяния заряженной частицы на кулоновском центре TD = v — = m V 8imG2mf In A va 47rng?g mJ In Л где liTitirnf,qtJqf - масса и заряд соответственно пробной частицы и частицы поля. In Л — Іп 53 - кулоновскии логарифм. Равенство массы и гравитационного "заряда" приводит к более простому выражению. Подставляя значения постоянных, получим следующую формулу V3 TD as 2 X 10е—2Млетп] rn2n где VCJ выражена в километрах в секунду, та в массах Солнца, an- число звезд в кубическом парсеке. Взяв га = 1, v = 100, и из средней плотности вещества в галактике 10 2С кг/м3 получаем оценку ті 0.27. Тогда время релаксации TD « 2 х 10илет
Для типичных параметров эллиптической галактики Т& 1014 лет, для дисковых галактик даже в центральной области время релаксации обычно превышает 1011 лет. Для сравнения, время жизни наиболее ярких звезд 106 — 10т лет, возраст галактики порядка 1010 лет.
Этот результат можно получить и другим путем. Для ограниченной одно-компонентной невращающейся системы радиуса R из N звезд скорость хаотического движения частиц определяется из теоремы вириала с2 = GNm/R, и время пролета тс = Rfc. Таким образом, получаем оценку [3] TD _ N тс InN отсюда видно, что большие звездные системы (в типичной галактике N 1010 -і-1012) должны рассматриваться как бссстолкновительные Аналогичные оценки справедливы и для пылевых частиц в протопланетном диске.
Таким образом, самогравитирующая система с достаточно большим количеством частиц может считаться бесстолкновителыгой. В такой системе, когда парных столкновений играют пренебрежимо малую роль, релаксация происходит за счет резонансов, важнейшим примером которых является затухание Ландау (применительно к гравитирующим системам это явление называется резонансом Линдблада).
Уравнение Пуассона
Как было отмечено в разделе 2.2,1, основную трудность в моделировании протопланетного диска представляет решение уравнения Пуассона. Перечисленные в разделе 1 особенности задачи накладывают следующие требования на метод решения уравнения Пуассона: решение в цилиндрической системе координат способность решать системы большой размерности за приемлемое время возможность эффективного распараллеливания Желательно, чтобы весь вычислительный эксперимент по моделированию протопланетного диска (около десяти тысяч временных шагов) выполнялся за время, не большее, чем 1 сутки. Это означает необходимость использования суперЭВМ и разработки параллельных методов решения уравнения Пуассона.
Гравитационный потенциал, в котором движутся частицы, можно разбить на две части: ф =$! +Ф2, где Фі - гравитационный потенциал массивной звезды, которая дает потенциал точечной массы . М0 Фі = -. г Здесь MQ - масса звезды. Вторая часть потенциала Фд определяется собственной массой пылевой среды и удовлетворяет уравнению Пуассона, которое в цилиндрических координатах имеет вид: 1 д2Ф2 д2Ф2 Здесь p(r, ifyz) - плотность пылевой среды, которая связана с функцией распределения /, входящей в уравнение Власова, соотношением р — / fdv
Граничные условия при г = Rm.ax и z = Zjnaj. задаются из точного решения уравнения Пуассона в случае, когда вся масса диска сосредоточена в точке (0,0,0). Точность задания этого условия определяется расстоянием от границы до диска.
По углу tp ставятся периодические граничные условия Ф2(гт27г,2) = 3( ,0,2). При г = 0 граничное условие имеет вид
При решении уравнений диск считается бесконечно тонким, в этом случае плотность вещества во всем объеме равна нулю (р = 0), все вещество находится в плоскости диска и рассматривается поверхностная плотность диска o(r,ip). В результате Z = 0 происходит разрыв нормальной производной потенциала, который дает граничное условие для определения потенциала Фз
Согласно [1] при наличии массивного центрального тела (масса центрального тела более чем в 5 раз превышает массу газового диска) движение газа может считаться двумерным. Таким образом, уравнения газовой динамики можно записать без учета производных по -координате Система уравнений газовой динамики в дивергентной форме имеет вид: + У(«пО=0 дії , 1 - — + (CV) и = —Vp + F ut О (4) { at є + ) + (tfV) (є 4- ) = —V (pfi) + + Fu- і«ДГ 2 у \ I j сг а а здесь Q - внешние источники тепла, к - коэффициент теплопроводности, F = УФ - внешние силы (сила гравитации и сила трения).
Давление р и внутренняя энергия є зависят от плотности и температуры: р = р( т,Т), є = є( г,Т), Уравнения состояния идеального газа: р = pRT, є = є(Т) Если є линейно зависит от Т, то (5) с = Ср RT 7-1,7 где R - универсальная газовая постоянная, Ср - теплоемкость газа при постоянном объеме, Си - теплоемкость газа при постоянном давлении.
Для интегрирования уравнения Власова-Лиувилля в ситуации, описанной в разделе 1 подходит только метод частиц в ячейках [96, 97\ 65, 9]. Согласно этому методу все вещество разбивается на модельные частицы, движение которых определяет эволюцию всей системы.
Основываясь на [9], приведем описание основной идеи метода- Пусть задача записана в абстрактной операторной форме; (б) Решение задачи (б) представляется в виде следующей интерполяционной формулы: N q = Sq,R(u,u,(t)). (7) J=i Этот переход называют разбиением среды на молельные частицы. Функция ii(u,v) называется ядром модельной частицы. Представление (7) позволяет свести решение задачи (б) к интегрированию динамической системы частиц.
Для протопланетного диска уравнение (6) - это уравнение Власова , и оператор А имеет вид A = v— —, (8) or т ov где г, г/, m - координата, скорость и масса частицы, Ф - гравитационный потенциал. Уравнения движения отдельных частиц имеют вид: dt u Эти уравнения характеристик для уравнения Власова. Решение уравнений движения для модельных частиц осуществляется по известной схеме Бориса [10]» которая состоит в решении уравнений движения в декартовых координатах с последующим пересчетом координат и скоростей частиц в цилиндрические координаты: Х = Тт + TVri У = TVp, rm+1 = л/х2-\-у2, sina = j//rm+l, cos a = x/r 1, ipm = pm+a. u+1 — vr cos a + 5 sin a, v+1 = i cosa — vrsina, VV » Здесь vTi Vp, T tp — скорости и координаты отдельных частиц, г - временной шаг, FT, F p - компоненты силы гравитации, действующей на частицу, интерполированные из узлов сетки в местоположение частиц.
Уравнения газовой динамики решаются методом крупных частиц Бслоцерковского-Давыдова [08]. Этот метод был выбран потому, что он хорошо согласуется с методом частиц в ячейках, как указано в [9]. Вместе с системой (4) решается уравнение температуры газа для того, чтобы избежать накопления ошибок при расчете полной энергии газа. Причина возникновения этих ошибок заключается в том, что в вакууме уравнения газовой динамики не могут быть записаны в виде (4). С вычислительной точки зрения трудность состоит в том, что уравнения содержат множитель о--1 {а - плотность газа), что приводит к численной неустойчивости при небольшом значении а.
Разностная схема для эйлерова этапа метрда крупных частиц
Время и точность решения уравнения Пуассона в значительной степени зависят от граничных условий. Поэтому необходимо разработать такой метод расчета граничных условий, который обеспечивал бы максимальной быстрое решение уравнения Пуассона при заданной точности. Расчетная область имеет вид кругового цилиндра: О г Rg 0 z Zg (59) По угловой координате существует условие периодичности Ф(0) = Ф(2тг) (GO)
В случае бесконечно тонкого диска правая часть уравнения (3) заменяется на 8{z)o{r,ip), где S(z) - дельта- функция Дирака, a{r, ip) - поверхностная плотность.
На поверхности диска существует разрыв нормальной производной потенциала дФ дФ — -— =2тг7а 61 OZ г-И-0 OZ г- -0 где о" - поверхностная плотность. Таким образом, на поверхности диска задается граничное условие второго рода. Для завершения постановки задачи необходимо еще задать потенциал на границе расчетной области.
Граничные условия для решения уравнения Пуассона естественным образом задаются на бесконечном расстоянии от диска: Ф(г) = О
Однако при численном решении уравнения Пуассона бесконечная расчетная область не может рассматриваться. Поэтому необходимо ограничивать расчетную область некоторым цилиндром, на границе которого задается граничное условие первого рода: Ф(г,2,)-Ф,(г) Ф(П9,г) = Ф2(2) здесь Rg и Zg - радиус и высота расчетной области.
При выборе размера расчетной области необходимо руководствоваться следующими соображениями: Граничные условия должны обеспечивать точное вычисление потенциала внутри расчетной области Алгоритм вычисления граничных условий не должен быть слишком трудоемким по сравнению с методом решения уравнения Пуассона
В соответствии с первым требованием границу расчетной области необходимо придвинуть как можно ближе к границе диска. Но в этом случае граничные условия необходимо будет рассчитывать заново на каждом временном шаге. Это приводит к увеличению времени счета. Таким образом, необходимо определить размер расчетной области исходя из соотношения времени расчета граничных условий и точности решения уравнения Пуассона.
В качестве первого варианта решения этой проблемы можно отодвинуть границу настолько далеко, что особенности распределения плотности на диске окажутся несущественными, В этом случае граничные значения можно вычислить по формуле потенциала точечного заряда: М , М Ф Ц " у/Щ ФіМ = --7-= (ад = здесь М - масса диска. Перерасчет потенциала при этом не нужен, так как масса диска не изменяется- Но у этого подхода есть два недостатка- Во-первых, требуется много дополнительных узлов сетки и, во-вторых, точность расчета для несимметричного диска достаточно мала.
Для того, чтобы учесть иесимметрига диска и сократить число дополнительных узлов сетки, граничные условия следует вычислять с помощью интеграла Ф(г0,р(ь-го) — - J rdr J —j-— py(rocos po - rcos p)2 + (rosinyjo -rsinip)2 4- o о о
Прямое его вычисление для каждой точки границы представляется слишком трудоемким. Если число узлов по каждому измерению равно JV, то сложность вычисления граничных условий будет Q(N4)
Более простой способ вычисления граничных условий - разложение в ряд Тейлора, был предложен В.Н.Снытниковым. М 1 Ф(г0, Vo, ) = -— - 2 5 № + » + - Щ здесь Ix, Iyj /І, Jo - осевые моменты инерции, В этом варианте граничные условия вычисляются относительно быстро {0(N2)), и появляется возможность придвинуть границу расчетной области за счет перевычисления граничных условий на каждом шаге.
Первый шаг - это дискретизация уравнения Пуассона конечно-разностным методом. Далее следует быстрое преобразование Фурье по угловой координате. Именно на этом этапе за счет линейной независимости Фурьс-гармоник удастся получить существенно параллельный метод решения. Это означает полное отсутствие межпроцессорных коммуникаций на стадии вычисления отдельных гармоник. В заключение раздела рассматриваются свойства двумерной СЛАУ, возникающей после прсобразооаиия Фурье, и описывающей отдельную гармонику. Эти свойства будут использованы при построении численных методов решения СЛАУ в разделе 4.3.3.
В цилиндрической системе координат вводится равномерная сетка. Дискретизация уравнения Пуассона на этой сетке дает систему линейных алгебраических уравнений
Декомпозиция области
Алгоритм, реализующий метод частиц в ячейках, обладает высокой степенью внутреннего параллелизма, так как частицы непосредственно друг с другом не взаимодействуют и для вычисления потенциала в слое сетки требуются значения только в двух соседних слоях. Тем самым выполнено условие линейности алгоритма [106], необходимое для его реализации на мультикомпыотерах. Несмотря на это, метод является сложным для эффективной параллельной реализации. Одним из главных вопросов - это правильное распределение частиц и массивов сеточных переменных между процессорными элементами (ПЭ).
Первая проблема состоит в минимизации объемов пересылок данных между ПЭ при реализации метода частиц на суперкомпьютерах с распределенной памятью. Например, при решении уравнения Пуассона необходимо знать значения потенциала в граничных узлах сетки, которые находятся в других процессорных элементах. Чтобы вычислить значение плотности в конкретном узле, необходимо знать координаты и массы всех частиц, находящихся в ячейках, одной из вершин которых является рассматриваемый узел. Однако частицы, принадлежащие соседним или даже одной и той же ячейке, могут храниться в разных процессорных элементах. Из этого следует, что пересылок данных избежать невозможно.
Другая проблема распараллеливания метода частиц состоит в трудности обеспечения равномерной загрузки процессорных элементов, поскольку первоначально равномерно распределенные в пространстве частицы через несколько шагов моделирования могут образовывать неравномерное распределение. Минимальный объем пересылок во время решения уравнения Пуассона и вычисления плотности достигается в том случае, когда описания частицы и значения потенциала в вершинах соответствующих ячеек хранятся в памяти одного процессорного элемента. В этом случае лучше всего, если частицы распределяются между ПЭ в зависимости от их размещения в пространстве моделирования, т.е. когда частицы, находящиеся в соседних ячейках, располагаются в одном и том же, либо в соседних ПЭ.
Первым способом решения этих проблем является дублирование в памяти каждого ПЭ сеточных значений потенциала во всем пространстве моделирования. Тогда множество частиц может быть распределено между ПЭ вне зависимости от юс расположения в пространстве моделирования и не потребуется передавать частицы из одного ПЭ в другой в ходе моделирования. Это так называемая Лагранжева декомпозиция. При этом существенно ограничивается объем решаемой задачи, так как массивы значений потенциала имеют большой объем (обычно это трехмерные массивы). Кроме того возникает необходимость передавать большие объемы данных между ПЭ, так как для расчета движения частиц в каждом ПЭ необходим потенциал вообще говоря, во всем пространстве моделирования. При этом частицы принадлежащие соседним или одной и той же ячейке могут храниться в разных ПЭ, поэтому на каждом пременном шаге необходимо суммировать значения средней плотности во всех узлах сетки, вычисленные в каждом ПЭ.
Второй способ распараллеливания метода частиц - Эйлерова декомпозиция - заключается в следующем: область моделирования разбивается на слои постоянного размера в каком-либо направлении, включая гравитационный потенциал и частицы в соответствующих точках. Число этих слоев определяется числом ПЭ. Направление разбиения выбирается таким образом, чтобы количество ячеек и количество частиц в каждом слое было примерно одинаковым. Б нашей задаче таким направлением является направление по угловой переменной. Тогда в начальный момент времени процессорные элементы имеют равномерную загрузку, В ходе расчета частицы, двигаясь в пространстве моделирования, могут переходить из подобласти хранящейся в одном ПЭ в подобласть, хранящуюся в другом ПЭ. Это приводит к неравномерному распределению нагрузки процессоров и возникает проблема балансировки загрузки [104, 105].
Если гармоники потенциала равномерно распределены по процессорам, и СЛАУ для каждой из гармоник решается методом БПВР, то распределение загрузки процессоров будет таким, как показано па рисунке 20. Загрузкой процессора в этом случае означает суммарное число итераций по всем гармоникам, находящимся на этом процессоре.
Как видно из рисунка 20, процессоры загружены неравномерно. Так как СЛАУ для 0-й гармоники обусловлена хуже всего, то процессор, содержащий эту гармонику (процессор номер 1), всегда будет перегружен по сравнению с другими. Таким образом, время решения уравнения Пуассона ограничено снизу временем решения СЛАУ для 0-й гармоники потенциала. Эта граница достигается в предельном случае, если число процессоров равно числу гармоник потенциала.
Как показано на рисунке 15т время расчета 0-м гармоники потенциала при использовании многосеточного метода значительно сокращается (для сетки 511 х 511 в 5.4 раза). При этом загрузка процессоров становится более равномерной. Дополнительно расчет первой гармоники выделяется на отдельный процессор. На рисунке 21 показано ускорение параллельной программы решения уравнения Пуассона. Сравниваются два варианта: равномерное распределение гармоник по процессорам с решением СЛАУ для всех гармоник методом БПВР и комбинированный вариант, включаюший в себя: Расчет нулевой гармоники многосеточным методом Выделение нулевой гармоники на отдельный процессор Расчет всех остальных гармоник методом БПВР и равномерное разделение их по процессорам
