CC BY
32
Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 2 (1), с. 201-207
УДК 524.52, 532.5
ПАРАЛЛЕЛЬНЫЙ КОД ДЛЯ ТРЕХМЕРНОГО МОДЕЛИРОВАНИЯ ПРОЦЕССОВ КОСМИЧЕСКОЙ ГАЗОДИНАМИКИ*
© 2011 г. М.А. Ерёмин, В.Н. Любимов
Волгоградский госуниверситет
ereminmikhail@gmail .com
Поступила в редакцию 24.09.2010
Разработан параллельный код для численного моделирования процессов космической газодинамики с учетом тепловых процессов. Представлены результаты трехмерного моделирования с высоким разрешением неупругих столкновений облаков межзвездного газа Н I при различных наборах начальных параметров. Проведено детальное исследование масштабируемости программы, в результате которого установлено, что эффективность распараллеливания составляет 70-90%, в зависимости от параметров задачи.
Ключевые слова: космическая газодинамика, трехмерное моделирование, параллельные вычисления, межзвездная среда, облака Н I, TVD-схемы.
Введение
Моделирование процессов космической газодинамики всегда являлось чрезвычайно ресурсоемкой сферой применения численных расчетов. Процессы, протекающие в астрофизических системах, являются трехмерными и, как правило, существенно нестационарными. Исследование подавляющего числа астрофизических задач может быть выполнено только численно в силу ограниченности аналитических методов. Отметим, что физические условия в астрофизических системах часто приводят к образованию ударных волн и сильной турбули-зации межзвездной среды.
Для адекватного моделирования астрофизических систем необходимо использовать трехмерные численные схемы и высокое пространственное разрешение, что делает возможным получение качественно новых результатов.
Одной из классических задач космической газодинамики является задача о столкновении облаков Н I в многофазной межзвездной среде [1]. Укажем на тот факт, что только при большом пространственном разрешении возможно наблюдать развитие гидродинамических и тепловых неустойчивостей, играющих решающую роль в разрушении холодных плотных облаков H I. С учетом трехмерности задачи это приводит к высоким требованиям к производительно-
* Статья рекомендована к печати программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2010» (http:// agora.guru.ru/pavt).
сти вычислительных систем и делает невозможным решение данной задачи без использования технологий параллельных вычислений.
Некоторые физические аспекты столкновений облаков в рамках достаточно грубых приближенных моделей впервые были рассмотрены в работах [2, 3]. Изучение столкновений облаков с использованием численного моделирования было проведено целым рядом авторов (см., например, работы [4-6]). Наиболее полное исследование столкновений облаков в двумерных численных моделях приведено в работах [7, 8], в том числе и с учетом магнитного поля. Тем не менее целый ряд вопросов остался без рассмотрения, особенно в трехмерном случае.
Целями нашей работы являются:
1. Разработка и реализация трехмерной параллельной ТУБ-схемы для моделирования процессов космической газовой динамики.
2. Применение разработанного кода для исследования неупругих столкновений облаков H I в многофазной межзвездной среде при различных начальных условиях, включая нелобовой характер столкновений и неидентичность сталкивающихся облаков.
3. Исследование масштабируемости параллельного кода и выявление основных причин издержек распараллеливания.
Постановка задачи
Процесс столкновения облаков атомарного водорода Н I, помещенных в теплую межзвездную среду, может быть описан уравнениями ньютоновской газовой динамики в одножид-
костном приближении. Предположим, что газ является идеальным и политропным с показателем адиабаты у = 5/3 .
В силу того, что масштаб Джинса намного больше размеров облаков X7 >> Яс, эффектами самогравитации газа пренебрегаем. Следует отметить, что на стадии компрессии облака испытывают существенное сжатие, что приводит к увеличению плотности, так что масштаб Джинса становится сравнимым с характерным вертикальным размером облака. Тем не менее в данной работе для упрощения модели самогравитация газа не рассматривается.
Система уравнений газовой динамики с учетом процессов нагрева, охлаждения и теплопроводности записывается в виде:
др
+ V • (ри) = 0,
+ V • (рм ® и) = ^р,
дг ^
дри
~дГ
дЕ ^ ^ ^
— + V • ([ Е + р]М) = -VI* (Т )VT ] -1( р, и).
дг
Здесь р,р,и = {и,V, м} - плотность, давление и вектор скорости газа соответственно,
(
Е = р
.2 Л
е + ■
и
2
- объемная энергия, к(Т) - ко-
мулах (2) учитывается через величины оа и ае, представляющие отношение классического теплового потока к насыщенному [15-17]:
ЯсЛа^гс _ каТ VT I
5р(квТ / то)312
Яс^с КТ5/2 ^Т|
\3/2
(3)
(4)
Яа 5р(квТ / то)
где т0 - масса протона, ка = 2.0 40—3
эрг/(е^см^К3/2) и ке = 5.6 -10-7 эрг/(схм^К7/2) -коэффициенты атомной и электронной теплопроводности соответственно.
Система уравнений газовой динамики (1) записывалась в безразмерной форме с использованием нескольких базисных величин. Перечислим наиболее важные размерные параметры задачи: Ь0 = 1 пк - характерный простран-
(1)
ственный масштаб, Т0 = 104 K - температура, и0 = 0.1 г/см3 - концентрация, сж0 «12 км/с -адиабатическая скорость звука, г0 «105 лет -
эффициент теплопроводности, Ь( р, и) =
= Л(Т)и2 — Г(Т)и - функция тепловых потерь, Л(Т) - функция объемного охлаждения, Г(Т) -функция объемного нагрева. Для идеального газа удельная энергия определяется выражением е = —р— . Связь между давлением и тем-
р(у—1)
пературой задается уравнением состояния: р = иквТ , где к в - постоянная Больцмана.
Функция охлаждения в наших численных моделях принималась отвечающей стандартному химическому составу межзвездной среды в окрестности Солнца (см., например, [9-11]). Расчеты показывают, что функция нагрева слабо зависит от температуры, поэтому мы полагали ее постоянной Г = 1.6 • 10—25 эрг/с.
Коэффициент теплопроводности учитывает как атомную, так и электронную теплопроводности, зависимость теплопроводности от температуры имеет следующий вид:
\кТ1/2 /(1 + а ), при Т < 104К,
*(Т) Ч 5/2/ 4 (2)
|кеТ5/2/(1 + сте), при Т > 104К.
При высоких температурах возникает эффект насыщения теплового потока, что в фор-
характерное время задачи.
Большое значение имеют равновесные параметры использованной физической модели. Укажем числовые значения следующих параметров, отвечающих тепловому равновесию: равновесная температура теплой фазы Те?1 = 8922 ^ равновесная температура холодной
фазы Те?2 = 62.5 К; концентрации и0 = 0.1 г/см3
соответствует равновесное давление межзвездной среды рщ /кв = 8922 К/см3.
Разработка и реализация кода
Численная схема
Для компьютерного моделирования неупругих столкновений облаков в межзвездной среде мы реализовали явную численную схему для (1) в декартовой системе координат. При построении численной схемы использовался метод расщепления по физическим процессам, согласно которому численное решение строится как решение уравнений, описывающих различные физические процессы.
Гидродинамические процессы, описываемые левой частью уравнений (1), представляют собой уравнения в частных производных гиперболического типа, для которых разработаны схемы неубывания полной вариации ТУБ [12], эффективно подавляющие нефизичные осцилляции численного решения на скачках. Реализованная пространственно не расщепленная
1Й1Ур|/р], 1=2.0
І8[|Ур[/р], 1=10.8
7
6
5
^ч 4 ш
ҐІ % \ к ^ \ 2
1 1
Рис. 1. Численный теневой шлирен (lg|Vр| / р ) для трехмерной модели лобового столкновения идентичных облаков (сечение плоскостью 1 = 0) в различные моменты времени
TVD-схема относится к типу М^^ [13], имеет третий порядок аппроксимации по пространству в областях гладкого течения и первый на скачках. За счет применения алгоритмов пересчета типа Рунге - Кутта реализованная численная схема обладает вторым порядком по времени.
Учет тепловых процессов, описываемых правой частью системы (1), сводится к коррекции тепловой энергии на каждом временном слое после расчета гидродинамических параметров течения. Изменение внутренней энергии за счет тепловых процессов описывается уравнением на внутреннюю энергию, которое решалось численно как с помощью метода подшагов, так и с помощью метода, предложенного в работе [7].
Распараллеливание
Для распараллеливания численного кода использовалась библиотека МР! [14]. Нами была выбрана линейная декомпозиция расчетной области на равные по объему участки. При таком делении расчетной области происходит уменьшение количества вызовов функций обмена данными. Однако при высоком разрешении и большом количестве процессов площадь границ между подобластями становится слишком большой и такой метод уже непрактичен. В нашем же случае - при использовании несколь-
ких десятков ядер - линейная декомпозиция является оптимальным вариантом.
Обмен границами производился в два этапа: сначала параллельно обменивались данными первый и второй, третий и четвертый и т. д. процессы, затем - второй и третий, четвертый и пятый и т. д. процессы. Таким образом, время обмена границами не зависит от количества подобластей. Реализованный способ обмена границами свидетельствует о правильности выбора линейной декомпозиции расчетной области, так как для трехмерной декомпозиции обмен границами пришлось бы проводить как минимум в шесть этапов (по количеству стыкующихся с соседними подобластями граней у каждой подобласти). Отметим также, что в программе используются асинхронные функции передачи данных MPI. Это позволяет во время обмена границами продолжать расчет в некоторой части подобласти, что приводит к минимизации задержек при обмене границами.
Результаты
Моделирование
Во всех численных моделях предполагалось, что в начальный момент времени температура и концентрация газа, отвечающие теплой фазе,
равны Тк = 104 К и пк = 0.1 г/см3, температура
облаков Н I Т = 80 К, концентрация газа в них в 100 раз больше, чем в межоблачной среде: П = 10 г/см3. Кроме того, в начальный момент
теплая и холодная фазы межзвездной среды находятся в состоянии механического равновесия, то есть рк = рс. Облака начинают движение навстречу друг другу с относительным числом Маха М = ис/с^ , где ск - адиабатическая скорость звука в межоблачной среде, ис - скорость движения облака.
Радиус сталкивающихся облаков варьировался от 0.5 пк до 1.5 пк. Нами были рассмотрены различные варианты столкновений: лобовые столкновения идентичных облаков, лобовые столкновения неидентичных облаков и нелобовые столкновения.
Расчеты проводились с разрешением до 12003 ячеек. Основной моделью является модель трехмерного лобового столкновения облаков Н I. На рис. 1 показаны основные стадии эволюции: а) стадия компрессии или сжатия, когда ударные волны, образовавшиеся при столкновении, распространяются наружу;
Ь) стадия перерасширения, на которой волны разрежения, распространяясь назад, формируют центральную область низкого давления и плотности; на этой стадии образуется выброс, подверженный действию неустойчивостей типа Кельвина-Гельмгольца; с) стадия коллапса, когда расширение останавливается давлением внешней окружающей среды; на этой стадии развивается неустойчивость Релея-Тейлора; d) стадия разрушения облаков.
Следует отметить, что для адекватного моделирования лобовых столкновений идентичных облаков необходимо использовать модели с высоким пространственным разрешением (50100 ячеек на радиус облака). Только в этом случае возможно образование филаментных структур, приводящих к фрагментации облаков. В противном случае происходит сильная деградация численного решения, вследствие чего в результате слияния двух сталкивающихся облаков возможно формирование нового облака, что как раз и наблюдалось в работе [7].
Также интерес представляет трехмерная модель нелобового столкновения облаков (рис. 2). Идентичные облака с радиусом 0.5 пк двигают-
Таблица 1
Эффективность программы
Количество MPI-процессов Время выполнения программы (сек) Время обмена границами (сек) Доля времени обмена границами (%) Эффективность программы (%)
1 12713 0 0 100
4 3285 8.4 0.26 96.76
20 733 21.1 2.88 86.72
40 470 45.85 9.76 67.63
Таблица 2
Зависимость времени счета от распределения MPI-процессов по узлам
Количество MPI-процессов Количество задействованных узлов кластера Время выполнения (сек)
8 1 4324
8 2 2538
8 4 2620
ся навстречу друг другу (Мх = М2 = 1.5) таким
образом, что центры облаков смещены на расстояние, равное одному радиусу. В результате столкновения образуются макроскопические вихревые филаменты с характерным размером несколько парсек.
Масштабируемость
Все расчеты и тесты проводились на соро-каядерном кластере, состоящем из пяти узлов, объединенных высокоскоростной шиной InfiniBand. Каждый узел содержит по два процессора Intel Xeon и 8 ГБ оперативной памяти.
Эффективность программы определялась как отношение ускорения программы S при использовании p процессоров к числу процессоров: E = Sp /p .
Ускорение программы вычисляется как отношение времени выполнения параллельной программы T к времени выполнения последовательной программы T1: Sp = Tp / T1.
Для измерения времени, затрачиваемого на обмен границами, в тестах использовалась версия программы с блокирующей пересылкой данных. Результаты тестов представлены в табл. 1.
Как видно из таблицы, время обмена границами сильно зависит от количества процессов. Данный факт требует объяснения, так как он вступает в противоречие с утверждениями пункта «Распараллеливание», в котором обосновывается выбор линейной декомпозиции расчетной области.
В тестовых расчетах задействовались все 5 узлов кластера и процессы равномерно распределя-
лись по всем узлам. В результате при большем количестве MPI-процессов количество задействованных ядер на узле возрастает, что приводит к увеличению конкуренции за общие ресурсы: оперативную память и канал передачи данных.
Для проверки зависимости скорости счета от распределения процессов по узлам мы провели тесты при одинаковом количестве MPI-процессов, но различном их распределении по узлам (табл. 2).
Во всех тестах количество используемых ядер равно восьми, но они располагаются на одном, двух или четырех узлах. Таким образом, при сохранении суммарной производительности вычислительной системы происходит изменение степени конкуренции за общие ресурсы. Время счета меняется до двух раз, причем самый медленный вариант наблюдается при расположении всех MPI-процессов на одном узле - когда максимальна конкуренция за общую оперативную память. Из этого следует вывод, что из общих ресурсов оперативная память является наиболее критичным источником замедления в нашем случае. Для уменьшения влияния данного фактора на оценку эффективности программы мы провели тесты, сохраняя число MPI-процессов на узле равным числу ядер (табл. 3).
Следует заметить, что при использовании от 16 до 40 процессов не наблюдается роста времени обмена границами, несмотря на то, что имеет место возрастание конкуренции процессов за канал передачи данных InfiniBand - десятков ядер явно недостаточно для проявления эффекта насыщения интерконнекта. Напротив, при переходе от 8 к 16 процессам (от использования одного к двум узлам кластера) имеет место значительный рост времени обмена грани-
Таблица 3
Эффективность программы при использовании всех ядер на каждом узле
Количество МР1-процессоров (количество задействованных узлов) Время выполнения программы (сек) Время обмена границами (сек) Доля времени обмена границами (%) Эффективность программы (%)
1 (1) 24913 0 0 100
8 (1) 4359 37.72 0.87 71.44
16 (2) 2095 112.21 5.36 74.32
24 (3) 1447 94.15 6.51 71.74
32 (4) 1162 123.27 10.61 67
40 (5) 921 93.67 10.17 67.63
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%,
-♦-Тест 1: эффективность -♦-Тест 1: издержки на обмен границами —Тест 2: эффективность "*“Тест 2: издержки на обмен границами
т К 1 і і і і і і і і і ■ И • _ її — ^
количество МРІ процессов
Рис. 3. Масштабируемость программы: тест 1 - равномерное распределение МР1-процессов по всем узлам кластера; тест 2 - по 8 процессов на каждом используемом узле
цами. Это вызвано появлением пересылок через канал передачи данных, помимо взаимодействий внутри одного узла.
На рис. 3 наглядно представлены результаты тестов в виде графиков. Во втором тесте в начале наблюдается резкое падение эффективности, связанное, как мы объяснили, с конкуренцией за оперативную память. Далее график принимает почти горизонтальный вид, и дальнейшее падение эффективности определяется только соотношением времени счета и времени обмена границами. Имеющееся в нашем распоряжении оборудование не позволяет провести тесты с более чем сорока ядрами и увидеть эффект насыщения реализованной параллельной программы.
Основные выводы
В заключение нашей работы сформулируем основные выводы:
1. Разработана и реализована трехмерная параллельная схема для численного моделирования процессов космической газодинамики.
2. Разработанный код обладает эффективностью около 70%.
3. Моделирование столкновений облаков Н I с использованием полученного кода показывает, что вне зависимости от характера столкновения происходит разрушение облаков с образованием филаментных структур.
Работа поддержана грантом ВолгУ 70-2011 -а/ВолГУ и грантами РФФИ № 11-02-97124-
р_поволжье_а, № 11-02-01332-а, а также ФЦП Рособразования (госконтракт П1248).
Список литературы
1. McKee C.F., Ostriker J.P. A theory of the interstellar medium - Three components regulated by supernova explosions in an inhomogeneous substrate // The Astrophysical Journal. 1977. Vol. 218. P. 148-169.
2. Stone M.E. Collisions between HI clouds. I. Onedimensional model // The Astrophysical Journal. 1970. Vol. 159. P. 277-292.
3. Stone M.E. Collisions between HI clouds. II. Two-dimensional model // The Astrophysical Journal. 1970. Vol. 159. P. 293-307.
4. Hausman M.A. Collisional mergers and fragmentation of interstellar clouds // The Astrophysical Journal. 1981. Vol. 245. P. 72-91.
5. Gilden D.L. Clump collisions in molecular clouds - Gravitational instability and coalescence // The Astrophysical Journal. 1984. Vol. 279. P. 335-349.
6. Lattanzio J.S, Monaghan J.J., Pongracic H., Scy-wartz M.P. Interstellar cloud collisions // Monthly Notes of the Royal Astronomical Society. 1985. Vol. 215. P. 125-147.
7. Miniati F., Jones T.W., Ferrara A., Ryu D. Hydrodynamics of cloud collisions in two dimensions: The fate of clouds in a multiphase medium // The Astrophys-ical Journal. 1997. Vol. 491. P. 216-232.
8. Miniati F., Ryu D., Ferrara A., Jones T.W. Mag-netohydrodynamics of cloud collisions in a multiphase interstellar medium // The Astrophysical Journal. 1999. Vol. 510. P. 726-746.
9. Spitzer L. Physics of fully ionized gases. New York: Interscience, 19б2.
10. Марочник Л.С., Сучков А.А. Галактика. М.: Наука, 1984.
11. Бочкарев Н.Г. Основы физики межзвездной среды. М.: МГУ, 1991.
12. Harten A. High resolution schemes for hyperbolic conservation laws // Journal of Computational Physics. 1983. Vol. 49. № 3. P. 357-393.
13. van Leer B. Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov's methods // Journal of Computational Physics. 1979. Vol. 32. № 1. P. 101-13б.
14. Антонов А.С. Параллельное программирование с использованием технологии MPI. М.: МГУ, 2004. 71 с.
15. Cowie L.L., McKee C.F. The evaporation of spherical clouds in a hot gas. I - Classical and saturated mass loss rates // The Astrophysical Journal. 1977. Vol. 211. P. 135-14б.
16. Balbus S. A., McKee C. F. The evaporation of spherical clouds in a hot gas. III - Suprathermal evaporation // The Astrophysical Journal. 1982. Vol. 252. P. 529-552.
17. Giuliani J. On the dynamics in evaporating cloud envelopes // The Astrophysical Journal. 1984. Vol. 277. P. б05-б14.
18. Ricotti M., Ferrara A., Miniati F. Energy dissipation in interstellar cloud collisions // The Astrophysical Journal. 1997. Vol. 485. P. 254-2б2.
19. Эндрюс Г.Р. Основы многопоточного, параллельного и распределенного программирования. М.: Издательский дом «Вильямс», 2003. 512 с.
PARALLEL CODE FOR THREE-DIMENSIONAL MODELLING OF COSMIC GAS DYNAMICS PROCESSES
M.A. Eremin, V.N. Lubimov
A parallel code has been developed for the numerical modelling of cosmic gas dynamics processes taking into account thermal processes. The results have been presented of high-resolution three-dimensional numerical modelling of inelastic interstellar HI cloud collisions with different sets of initial parameters. A detailed study of program scalability has shown that parallelization efficiency amounts to 70-90 % depending on the problem parameters.
Keywords: cosmic gas dynamics, three-dimensional modelling, parallel computing, interstellar medium, HI clouds, TVD schemes.
