CC BY
43
Вычислительные технологии
Том 12, № 3, 2007
ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ НА СУПЕРЭВМ МОДЕЛИ ГАЗОВОЙ КОМПОНЕНТЫ САМОГРАВИТИРУЮЩЕГО ПРОТОПЛАНЕТНОГО ДИСКА*
В. А. Вшивков, Г. Г. Лазарева, С. Е. Киреев, И.М. Куликов Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: lazareva@ssd.sscc.ru
A parallel realization of the numerical algorithms for the 3D simulation of transient processes in gravitating gas systems with self-consistent field in Cartesian coordinates is considered. The evolution of rotating gas is described by a system of gas dynamics equations and the Poisson equation for the gravitational potential. This paper contains a detailed analysis of the contribution of each part of the algorithm into the total calculation time. The main characteristics of the parallel algorithm are presented.
Введение
Эволюция самогравитирующих газовых систем вызывает большой интерес у астрофизиков [1]. Для изучения зарождения звезд и планет необходимо численное моделирование нестационарной и пространственно трехмерной динамики газа [2, 3]. В работе представлена параллельная реализация численной модели для трехмерного моделирования нестационарных процессов в гравитирующих системах с самосогласованным полем. Она основана на решении уравнения Пуассона для гравитационного поля и газодинамических уравнений. Газодинамическая часть представлена уравнениями для плотности, вектора скорости и энергии в трехмерной декартовой системе координат. Для последующего добавления расчета химических реакций система уравнений газовой динамики дополнена уравнением для температуры, полученным как следствие закона сохранения внутренней энергии. Трехмерность модели и нестационарность задачи выдвигают строгие требования к экономичности методов ее решения как в плане использования ресурсов вычислительной системы, так и в плане некритичного ограничения на отношение шага по времени и пространству. В последнее время благодаря бурному развитию вычислительной техники стали возможны
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00665), программы СО РАН по суперЭВМ, программы Президиума РАН "Происхождение и эволюция биосферы" № 18-2 (2006), программы Президиума РАН П-04 "Происхождение и эволюция звезд и галактик", программы Рособразования "Развитие научного потенциала ВШ" (проекты РНП.2.2.1.1.3653 и РНП.2.2.1.1.1969).
( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
ресурсоемкие расчеты с получением физически оправданных результатов для трехмерных программ. Суперкомпьютеры позволяют использовать большие объемы данных, на порядки повышать производительность вычислений, а как следствие, и точность. Появление параллельной версии программы для трехмерного моделирования динамики газовой компоненты протопланетного диска в гравитационном поле позволит создать условия для протекания химических реакций и получения органических и неорганических соединений.
Одна из проблем, с которой пришлось столкнуться в последние годы, связана с трудностями, возникающими при использовании многопроцессорных суперЭВМ. Эта проблема в первую очередь определяется сложностью адаптации алгоритмов решения задач, описываемых уравнениями газовой динамики [4-9] и уравнением Пуассона [10-15] на архитектуру многопроцессорных вычислительных систем с распределенной памятью. Вопросы, связанные с эффективным использованием высокопроизводительной вычислительной техники, носят системный характер и их решение невозможно без привлечения фундаментальной науки [16]. Заметим, что, используя такие подходы к параллельной численной реализации математических моделей газодинамики, можно проанализировать правомерность применения самих моделей, определить диапазон допустимых параметров и оценить точность получаемого решения. Применяя готовые пакеты программ [17, 18], невозможно провести такое исследование методов решения.
В рамках данной задачи разработаны параллельный алгоритм и параллельная версия программы для моделирования вращающегося газового облака в самосогласованном гравитационном поле с центральным телом, с учетом температуры газа, в трехмерной постановке в декартовых координатах. Особенностью задачи является неустойчивость по Адамару исходной системы уравнений (отсутствие непрерывной зависимости решения от начальных параметров задачи) [19]. Физическая постановка задачи требует допущения о правомерности использования системы уравнений газовой динамики в области с нулевой плотностью. Вследствие такого допущения возникают проблемы численной реализации этой системы на границе сред. Для описания вращающейся газовой компоненты диска в вакууме нужен особый подход к моделированию границы газ—вакуум, исключающий счетные осцилляции или размазывания границы. Для решения системы уравнений газовой динамики используется модифицированный метод крупных частиц Белоцерковского— Давыдова [20]. С целью контроля правильности и регуляризации решения в систему введено независимое уравнение для температуры, реализованное с помощью схемы стабилизирующей поправки. Решение уравнения Пуассона находится методом быстрого преобразования Фурье. В результате комбинации схем решения с различными свойствами построена численная модель, реализованная в виде программного комплекса для расчетов на суперЭВМ. Работоспособность созданной параллельной версии кода демонстрируется на решении задачи о динамике газовой компоненты протопланетного диска в гравитационном поле.
1. Постановка задачи
Рассмотрим систему уравнений газовой динамики, дополненную уравнением Пуассона для гравитационного потенциала:
др
— + ¿1у(ру) = 0,
dt
dpE
~дГ
pv + div(vpV) = —grad p — p grad Ф,
+ div(pEv') = —div(pv') — (p grad Ф, v) + div(x grad T), p = pRT, АФ = 4np,
где p — давление, p — плотность, v — скорость, E — полная энергия, Ф — гравитационный потенциал, х — коэффициент теплопроводности, R — универсальная газовая постоянная.
Для обезразмеривания переменных в качестве основных характерных параметров выбраны:
— расстояние от Земли до Солнца R0 = 1.5 • 1011 м;
— масса Солнца M0 = 2 • 1030 кг;
— гравитационная постоянная G = 6.67 • 10-11 Н-м2/кг2.
Соответствующие характерные величины скорости частиц v0, t0, потенциала Ф0 и поверхностной плотности а0 могут быть записаны как
GM0 30 / . R0 „ 1/6 ф 2 GM0 M0
—— = 30 км/с, ¿0 = — ~ 1/6 года, Ф0 = v0 = ^0 = •
R0 v0 R0 R0
Так как исходная задача неустойчива и ее постановка некорректна, вопрос о контроле правильности решения численной реализации задачи стоит особенно остро. В первую очередь важна уверенность как в правильности решения системы уравнений газовой динамики без учета влияния гравитационного поля, так и в правильности решения уравнения Пуассона. В качестве тестовой задачи для системы уравнений газовой динамики использовались тесты Годунова [21, 22] и классическая задача о разлете газового шара в вакуум. Решение уравнения Пуассона протестировано на аналитических решениях. Для контроля правильности решения используются законы сохранения массы, импульса, полной и внутренней энергии, а также теорема вириала [1].
На границах расчетной области для газодинамических величин поставлены однородные краевые условия второго рода, для гравитационного потенциала граничные условия задаются с помощью мультипольного разложения с использованием инерциальных моментов системы [23].
В расчетной области в начальный момент времени задается нулевое распределение плотности и скоростей. В центре области задается область эллипсоидальной формы, заполненная газом с заданным распределением плотности, температуры, скоростей и гравитационного потенциала. С течением времени газовое облако в вакууме разлетается, кол-лапсирует или остается неподвижным в зависимости от задания начальных данных.
2. Метод решения
Задача решена в трехмерной постановке в декартовых координатах. В расчетной области задана равномерная прямоугольная сетка размером Кх х Ку х Kz. В узлах сетки определены значения трех компонент скорости, в центрах ячеек — остальные параметры задачи. Всего при решении уравнений газовой динамики используются 18 трехмерных массивов размерности Кх х Ку х Kz: 11 массивов для хранения параметров газа (плотность, давление, температура, энергия, гравитационный потенциал, скорость и импульс) и семь массивов для хранения промежуточных данных между различными этапами вычислений.
За основу метода решения системы уравнений газовой динамики выбран метод Бело-церковского—Давыдова. Этот метод обеспечивает автоматическое выполнение законов сохранения массы, импульса и полной энергии, но является методом первого порядка точности. При выборе метода было принято во внимание, что схема первого порядка точности обладает схемной вязкостью, которая подавляет счетную дисперсию. Появление коротковолновых возмущений масштаба расчетной ячейки, численных флуктуаций плотности, обусловленных счетной дисперсией, представляется менее нежелательным, чем некоторое размазывание скачков функций. Отсутствие счетных флуктуаций плотности необходимо для расчета физических неустойчивостей гравитационного типа. Кроме того, очень сложна реализация схем более высокого порядка точности на границах с вакуумом. Метод Белоцерковского—Давыдова расщепляется на два этапа: эйлеров этап (рис. 1, п. 3 схемы), на котором отбрасываются члены конвективного переноса, и лагранжев этап (рис. 1, п. 5 схемы), на котором происходит конвективный перенос газодинамических величин. Для реализации первого этапа используется конечно-разностная схема с операторным походом [24]. Разностная схема, построенная на основе операторного подхода, характеризуется различием в способах аппроксимации производных, описывающих разные дифференциальные операторы. В трехмерном случае схема реализуется на компактном 27-точечном шаблоне. Особенность такой схемы состоит в выполнении множества операций осреднения. Как известно, метод Белоцерковского—Давыдова слабоустойчив. Но в силу неустойчивости физического решения задачи такое свойство метода не имеет решающего значения.
Уравнение для полной энергии можно записать в следующем виде:
Вычисление температуры необходимо для последующего расширения задачи — включения в численную модель блока химических реакций. Решение уравнения для температуры независимо от решения системы уравнений газовой динамики позволяет проводить регуляризацию решения с сохранением консервативности схемы. Для решения уравнения температуры (рис. 1, п. 2 схемы) использовалась схема стабилизирующей поправки [25], обладающая абсолютной устойчивостью. Суть метода состоит в расщеплении исходного уравнения на три, каждое из которых реализуется скалярной прогонкой по каждому из координатных направлений. Каждое из трех уравнений реализовано с использованием операторного подхода и решается на компактном 27-точечном шаблоне. Тестовые расчеты с использованием последовательного алгоритма показали [23], что регуляризация решения с помощью уравнения температуры позволяет получать решения на границе газ—вакуум.
Значения температуры вычисляются независимо с целью контроля выполнения законов сохранения как полной, так и внутренней энергии (энтальпии). Контроль осуществляется перенормировкой (рис. 1, п. 4 схемы) схемных скоростей переноса массы, импульса и удельной энергии на лагранжевом этапе метода Белоцерковского—Давыдова введением множителя у/2(Е — Т) /1. Такая перенормировка сохраняет направление вектора скорости, корректируя его длину. В случае, когда неизбежное для каждого численного метода накопление погрешности вычислений приводит к отрицательным значениям внутренней энергии, вычисляемые значения полной энергии корректируются таким образом, чтобы энтальпия стала нулевой. В результате корректировки численная модель обеспечивает решение на границе с областью нулевой плотности, в которой система уравнений газовой динамики неправомерна.
дТ
— + ^га^, V) + (7 — 1)Т
1
¿1у(хёга^) = 0.
р
1. Распределение значений плотности, импульса, удельной полной энергии, давления, температуры и гравитационного потенциала (М = 1,..., N — 1, N — полное число процессоров)
0 М N — 1
2. Решение уравнения температуры
Первый этап решения уравнения температуры: • скалярная прогонка в направлении оси ОХ с передачей граничных значений между процессорами Прямой ход
М — 1 М М + 1
Обратный ход
М — 1 М М + 1
Второй и третий этапы решения уравнения температуры:
• скалярные прогонки в направлении осей ОУ и OZ выполняются независимо без межпроцессорных обменов;
• обмен граничными значениями температуры
М—1
М
М + 1
3. Реализация эйлерова этапа метода крупных частиц:
• независимое вычисление импульса и удельной полной энергии;
• обмен граничного значения
4. Независимая коррекция значений схемных скоростей
5. Реализация лагранжева этапа:
• независимый перенос значений плотности, импульса и удельной полной энергии;
• сборка граничных значений
М—1
М
М + 1
6. Независимое вычисление скорости и давления
7. Решение уравнения Пуассона:
• перенос значений плотности из малой сетки в большую (обмены);
• вычисление граничных условий (сборка значений со всех процессоров);
• прямое БПФ по трем направлениям (межпроцессорные обмены);
• независимое вычисление по схеме;
• обратное БПФ по трем направлениям (межпроцессорные обмены);
• перенос значений гравитационного потенциала из большой сетки в малую (межпроцессорные обмены)
Рис. 1. Схема работы параллельного алгоритма на одном временном шаге
Результатом такой комбинации схем решения с различными свойствами является численная модель, которая обладает всеми необходимыми качествами. Использование слабоустойчивого метода Белоцерковского—Давыдова дает основание говорить о том, что физическая неустойчивость задачи не нарушается. Кроме того, метод хорошо описывает границу газ—вакуум. Введением абсолютно устойчивой схемы для реализации уравнения температуры удалось повысить устойчивость всей численной модели в целом.
Для получения распределения гравитационного потенциала в области моделирования на каждом временном шаге необходимо решать уравнение Пуассона. Дискретизация вторых производных в трехмерном уравнении Пуассона производится по 27-точечной схеме,
имеющем второй порядок аппроксимации по пространству: 1 38 4
^ (-уФгД! + д(Фг—1,к,1 + Фг+1,к,1 + Фг,к—1,1 + Фг,к+1,1 + Фг,к,1—1 + +
— 1 ,к—1,1 — 1,к+1,1 + Фг+1,к—1,1 + Фг+1,к+1,1 +
— 1,к,1—1 — 1,к,1+1 + Фг+1,к,1—1 + Фг+1,к,1+1 +
+ Ф г,к—1,1—1 + Ф г,к—1,1+1 + Ф г,к+1,1—1 + Ф г,к+1,1+1 ) +
+ 71(Фг— 1,к—1,1—1 + Фг—1,к—1,1+1 + Фг—1,к+1,1—1 + Фг—1,к+1,1+1 +
36
+ Фг+1,к—1,1—1 + Фг+1,к—1,1+1 + Фг+1,к+1,1—1 + Фг+1,к+1,1+1)) = Рг,к,1 • Значения потенциала на границе вычисляются по приближенной формуле:
М М
Ф(ж,у,г) = — + — (1х + /у + 4 - 31о).
г г 3
Формула получена с помощью мультипольного разложения с использованием моментов инерции системы:
I
0
(x2Ix + y2Iy + Z2 Iz) - + xzlxz + yzIyz)
4 = + y2)mfc' Jy= + zk)mfc> Iz = + > fc fc fc
Ixy = ^ xfcyfcmfc, Ixz = ^ xfczfcmfc, Iyz = ^ yfczfcmfc. fc fc fc
Здесь r = \Jx2 + y2 + z2, где x, y, z — координаты точки пространства относительно центра масс. В определении осевых и центробежных моментов инерции m^ — масса и Xfc, yfc, Zfc — координаты точки. Точное вычисление краевых условий для уравнения Пуассона суммированием потенциала по всей сетке привело бы к слишком большим временным затратам. Оценим максимальную относительную погрешность потенциала, полученного по приближенной формуле. Если масса системы сосредоточена не далее чем на расстоянии Ä! от центра масс, то легко показать, что максимальная относительная погрешность потенциала е в точке на расстоянии Я2 от центра масс будет равна
е = R1
я*'
Таким образом, чем дальше отодвинуть границу области при решении уравнения Пуассона, тем точнее будет результат. В связи с этим уравнение Пуассона решается в области большего размера, чем та, в которой движется газ. За пределами области движения газа плотность берется нулевой. Заметим, что для более простой формулы вычисления потенциала Ф = -M/r оценка погрешности будет выглядеть как е = Я2/Я2, что потребует для достижения той же точности отодвинуть границу значительно дальше, т. е. использовать значительно большие сетки.
Уравнение Пуассона решается методом преобразования Фурье. Этот метод имеет число операций порядка O(n log n) и является наиболее оптимальным в случае использования
равномерной пространственной сетки, так как очень прост в реализации и достаточно эффективен. Кроме того, данный подход требует использования всего одного массива как для значений плотности на входе, так и для значений потенциала на выходе решателя, что существенно при решении больших задач.
Время счета зависит в первую очередь от параметров вычислительной модели, т. е. от размера сетки. Временной шаг т выбирается постоянным. Тем самым исключается косвенная зависимость от физических параметров исследуемой системы (размера, массы, распределения плотности).
3. Выбор метода декомпозиции
Алгоритм, реализующий метод крупных частиц, обладает высокой степенью внутреннего параллелизма, так как имеет место только локальное взаимодействие соседних элементов и для вычисления потенциала в слое сетки требуются значения лишь в двух соседних слоях. Тем самым выполнено условие линейности алгоритма [26], необходимое для его реализации на мультикомпьютерах. Несмотря на это, метод является сложным для эффективной параллельной реализации. Один из главных вопросов — это правильное распределение массивов сеточных переменных между процессорными элементами. При распараллеливании решения уравнения температуры возникает проблема реализации скалярной прогонки, которая является существенно последовательным алгоритмом. Несмотря на наличие известных параллельных алгоритмов для реализации прогонки их использование дает только небольшой выигрыш в производительности, что с лихвой "компенсируется" сложностью реализации [27]. Для повышения точности решения уравнения Пуассона требуется отодвигать границу, а значит, использовать на этом этапе вычисления большей сетки. Необходимость переноса значений между сетками разного размера порождает дополнительные проблемы при распараллеливании.
Для распараллеливания алгоритма рассматриваемой задачи применен метод геометрической декомпозиции. Расчетная область размера Кх х Ку х Kz разрезается вдоль оси X (рис. 2, а) на слои размера (Кх/Р) х Ку х Kz, где Р — количество процессоров вычислительной системы. Все массивы "разрезаются" на подмассивы в соответствии с разрезанием области моделирования. Декомпозиция области осуществляется с перекрытием граничных точек соседних подобластей (рис. 2, б), поэтому размерность этих массивов равна (Кх/Р + 2) х Ку х Kz.
Рис. 2. Схема декомпозиции области моделирования (а), схема перекрытия граничных подобластей (б)
Такой способ декомпозиции области возможен, так как используемая модель динамики газовой компоненты протопланетного диска не разделяет область, занятую газовым шаром, и вакуум, т.е. позволяет считать среду однородной. Расчет уравнений газовой динамики ведется по схеме крупных частиц, являющейся схемой сквозного счета (без выделения особенностей), на компактном 27-точечном шаблоне. Для расчета гравитационного потенциала как в газовом облаке, так и в вакууме используется одно и то же уравнение Пуассона, в правую часть которого входит плотность газовой фазы, имеющая нулевое значение в вакууме.
Так как данные разрезаны на слои и обмены данными во время счета производятся только между соседними слоями, логическая топология вычислительной системы "линейка" достаточна для решения задачи. Слои декомпозированной расчетной области распределяются по компьютерам последовательно в соответствии с номерами компьютеров. Слой с наименьшими координатами узлов сетки располагается в нулевом компьютере, с большими координатами сетки — в первом компьютере и т. д., с наибольшими координатами сетки — в последнем компьютере.
Как сказано выше, уравнение Пуассона решается в области большего размера, чем расчетная область для решения системы уравнений газовой динамики. В дальнейшем область (и сетку), в которой решается уравнение Пуассона, будем называть "большой", а область (и сетку), в которой решаются уравнения газовой динамики, — "малой". На рис. 3, а малая область показана сплошной линией, а большая — пунктирной.
В параллельной реализации массив разрезается перпендикулярно оси ОХ на число частей, равное числу процессоров. Перекрытие областей при этом не требуется. Трехмерное параллельное быстрое преобразование Фурье выполняется с помощью процедуры из свободно распространяемой библиотеки ЕЕТШ [28]. Способ разрезания массива также задается библиотекой.
Сетка для решения уравнений газовой динамики вложена в сетку, на которой происходит решение уравнения Пуассона. На рис. 3, б изображено взаимное расположение большой и малой сеток в пространстве моделирования (проекция на плоскость XV), а также пример их распределения между процессорами. Способы разрезания этих двух сеток совпадают, но различие их размеров приводит к тому, что сеточные значения по-разному распределены между процессорами. Выбранный способ решения уравнения Пуассона не позволяет проводить вычисления отдельно на малой сетке. Поэтому до и после решения уравнения Пуассона необходимо выполнять перераспределение сеточных значений. Структура обменов при этом может существенно различаться в зависимости от размеров сеток
Рис. 3. Схема расположения вложенных сеток (а), схема декомпозиции большой и малой сеток (б), схема обменов при перераспределении сеточных значений (в)
и числа процессоров. На рис. 3, в показан пример схемы обменов для четырех процессоров. Числами 0, 1, 2 и 3 обозначены номера процессоров, на которых находятся соответствующие фрагменты сеток. Стрелками показаны передачи данных между сетками. Видно, что для переноса значений из одной сетки в другую требуются обмены между 0-м и 1-м, а также между 2-м и 3-м процессорами. Остальные блоки данных могут быть просто скопированы в памяти внутри 1-го и 2-го процессоров. Перераспределение сеточных значений — узкое место в представленном алгоритме, так как с ростом числа процессоров число обменов будет возрастать.
4. Ускорение и эффективность параллельного алгоритма
При разработке параллельного алгоритма важно знать потенциальные возможности ускорения вычислений и накладные расходы, связанные с организацией взаимодействий параллельных ветвей. Кроме того, важно знать показатели эффективности работы параллельного алгоритма на многопроцессорной системе, позволяющие сравнивать его с другими параллельными алгоритмами. Так как для получения значимого результата задачи динамики газовой компоненты протопланетного диска требуются большие объемы данных, важно оценить предел роста размеров задачи, после которого ускорение и эффективность параллельного алгоритма перестанут расти. Для обработки объема данных больше предельного нужно будет искать новые подходы к распараллеливанию алгоритма либо к решению всей задачи в целом.
Для оценки параллельного алгоритма наряду с временем вычислений будем рассматривать ускорение параллельного алгоритма (рис. 4, б)
иР = Т1
р Тр
и показатель эффективности распараллеливания на системе из Р компьютеров (рис. 4, в)
= Цр ■ 100 %.
Здесь ТР — время счета параллельного алгоритма на системе из Р компьютеров, Т1 — время счета параллельного алгоритма на одном процессоре. В расчетах использовались разрешения сетки 196 х 196 х 196 узлов для решения уравнения газовой динамики и 392 х 392 х 392 узла для решения уравнения Пуассона. Измерялось время, необходимое для расчета десяти временных шагов (рис. 4, а).
На рис. 5, а видно, что наибольшая вычислительная нагрузка приходится на решение уравнений газовой динамики. Подробно проанализируем вклад в общее время счета того времени, которое требуется многопроцессорной системе для численного решения уравнения температуры, а также эйлерова и лагранжева этапов метода крупных частиц при решении системы уравнений газовой динамики. Как видно из диаграммы зависимости времени счета от размера вычислительной системы (рис. 5, б), при малом числе процессоров наибольшую долю времени занимает расчет импульса и удельной энергии на эйлеровом этапе (см. рис. 1, п. 3 схемы). Операторный подход [24], примененный для численной реализации первого этапа схемы Белоцерковского—Давыдова, несомненно, требует для своей
реализации большого количества обращений к функциям вычисления разностных аналогов дифференциальных операторов. Это является существенным недостатком выбранного метода, поэтому возникает необходимость поиска более экономичных способов получения инвариантных относительно поворота решений задач гравитационной физики.
Следующий по количеству затрачиваемого времени — этап решения уравнения температуры. Метод прогонки, использующийся на данном этапе, является существенно последовательным алгоритмом. В результате при расчете на большом числе процессоров этот этап превосходит по времени счета все другие этапы.
На рис. 5, в представлены времена счета различных этапов решения уравнения Пуассона в зависимости от числа процессоров. Из графиков видно, что с ростом числа процессоров все этапы этой части алгоритма, за исключением операции перераспределения данных между сетками, дают ускорение.
В приведенном выше исследовании во всех случаях использовалась задача одного и того же размера. Следовательно, с ростом числа процессоров нагрузка на каждый отдельный процессор и подсистему памяти узла сокращалась. Полезным может оказаться исследование эффективности алгоритма на задачах, размер которых обеспечивает примерно одинаковую нагрузку на каждый процессор. На рис. 6 представлены результаты
Т
1600-
и
и-
16 32
1 2 4
а б в
Рис. 4. Время работы (а), ускорение (б) и эффективность (в) параллельного алгоритма
Рис. 5. Относительное время работы различных этапов алгоритма: а — время выполнения всей задачи: 1 — решение уравнений газовой динамики, 2 — решение уравнения Пуассона, 3 — другие вычисления; б — время решения уравнений газовой динамики: 1 — решение уравнения температуры, 2 — эйлеров этап, 3 — лагранжев этап, 4 — другие вычисления; в — время решения уравнения Пуассона: 1 — вычисление граничных значений, 2 — перераспределение данных между сетками, 3 — выполнение процедуры БПФ, 4 — локальные вычисления по схеме, 5 — другие вычисления
а
в
такого тестирования, когда на каждый процессор приходится около 750 МБ данных. Размер сеток при этом варьировался от 156 х 156 х 156 для малой и 312 х 312 х 312 для большой при запуске на одном процессоре и до 480 х 480 х 480 для малой и 960 х 960 х 960 для большой при запуске на 32 процессорах. Времена счета представлены на рис. 6, а. Ускорение на рис. 6, б вычислялось по следующей формуле:
и' = ТР
и = Тр '
где Тр — время счета задачи на Р процессорах. Величина Т'Р — это гипотетическое время счета на одном процессоре задачи, размер которой равен размеру задачи, посчитанной на Р процессорах. Эффективность на рис. 6, в вычислена по обычной формуле:
^ = -Р ■ 100 %.
Сравнивая рис. 4 и 6, можно увидеть, что с увеличением объема задачи на малом числе процессоров эффективность несколько ухудшилась. Но уже на 32 процессорах эффективность распараллеливания на большой задаче оказалась более высокой, чем на задаче
абв
Рис. 6. Время работы (а), ускорение (б) и эффективность (в) параллельного алгоритма при одинаковой загрузке процессоров
а
в
Рис. 7. Относительное время работы различных этапов алгоритма при одинаковой загрузке процессоров: а — время выполнения всей задачи: 1 — решение уравнений газовой динамики, 2 — решение уравнения Пуассона, 3 — другие вычисления; б — время решения уравнений газовой динамики: 1 — решение уравнения температуры, 2 — эйлеров этап, 3 — лагранжев этап, 4 — другие вычисления; в — время решения уравнения Пуассона: 1 — вычисление граничных значений, 2 — перераспределение данных между сетками, 3 — выполнение процедуры БПФ, 4 — локальные вычисления по схеме, 5 — другие вычисления
небольшого объема. То есть на большой задаче по сравнению с малой доля независимых вычислений оказалась больше. Рисунок 7 показывает, какие изменения во времени вычислений претерпели различные этапы алгоритма.
Видно, что, как и было отмечено ранее, узким местом при решении уравнений газовой динамики является решение уравнения температуры методом скалярной прогонки по трем направлениям. На этапе решения уравнения Пуассона увеличение размера задачи и числа процессоров хуже всего сказывается на операции быстрого преобразования Фурье. Значительное замедление при этом связано с выполнением операции поворота трехмерного распределенного массива относительно оси разрезания. Немонотонный характер графика объясняется различием в реализациях библиотечной операции БПФ для различных размеров массива. Кроме того, при увеличении числа процессоров растет время перераспределения данных между большой и малой сетками.
5. Результаты применения параллельного алгоритма
Результаты численного моделирования динамики газовой компоненты протопланетного диска в гравитационном поле в трехмерном случае с использованием описанного выше алгоритма приведены на рис. 8.
В начальный момент времени газ равномерно распределен в области, ограниченной шаром радиуса Я = 2.0. Начальная угловая скорость задается формулой:
Уф = ^(1 + сое 3ф).
Остальные компоненты скоростей равны нулю.
В ходе расчета газовое облако приобретает конфигурацию в виде трех рукавов плотности. Расчет на сетке 64 х 64 х 64 (128 х 128 х 128 — для решения уравнения Пуассона)
Рис. 8. Проекции распределения плотности газа на экваториальную (а и б) и меридиональную (в и г) плоскости при использовании сетки 64 х 64 х 64 (а и в) и 512 х 512 х 512 (б и г) в момент времени Ь = 2.4
проводился на одном процессоре. Полученное распределение плотности хорошо передает качественный характер решения. Многопроцессорная вычислительная система позволяет получить распределение плотности в газовом облаке с детальным описанием областей ее высоких градиентов. Расчет на сетке 512 х 512 х 512 (800 х 800 х 800 — для решения уравнения Пуассона) проводился на 32 процессорах. Используемый объем данных во втором случае составил более 21 Гбайт. Результаты трехмерных расчетов подтвердили эффективность разработанного численного метода решения задачи.
6. Выводы
В статье рассмотрены проблемы разработки параллельной программы для моделирования динамики газовых структур на многопроцессорных компьютерах с распределенной памятью. Созданная параллельная реализация дает адекватные результаты для трехмерной модели гравитационной газовой динамики. Одна из главных трудностей при этом состоит в вычислении гравитационного потенциала, который задается трехмерным уравнением Пуассона. При решении уравнения Пуассона необходимо отодвигать границу, что приводит к необходимости использования большей сетки. Эта сетка при разрезании является не согласованной с другими сетками, поэтому необходимо выполнять перераспределение сеточных значений между различными этапами алгоритма. Скалярная прогонка, используемая при численной реализации уравнения температуры, является существенно последовательным алгоритмом. Для ее эффективной параллельной реализации вдоль оси ОХ потребовалось использовать дополнительный массив. В результате на данном этапе алгоритма при использовании нескольких процессоров также было получено ускорение. Тестовые вычисления на суперкомпьютере МВС-1000М показали, что при решении системы уравнений газовой динамики линейного ускорения не происходит, так как имеет место насыщение в результате обменов. Решение системы уравнений газовой динамики совместно с уравнением Пуассона дало 12-кратное ускорение на 32 процессорах.
Реализация трехмерной модели гравитационной газовой динамики на многопроцессорной вычислительной системе позволило использовать для расчетов подробные сетки (512 х 512 х 512), а следовательно, дало возможность решать более широкий класс задач.
Список литературы
[1] Bertin G. Dynamics of Galaxies. Cambridge Univ. Press, 2000.
[2] Снытников В.Н., Плрмон В.Н. Жизнь создает планеты? // Наука из первых рук. 2004. № 0. C. 20-31.
[3] Снытников В.Н., Вшивков В.А., ДудниковА Г.И., Никитин С.А., Пармон В.Н., Снытников А.В. Численное моделирование гравитационных систем многих тел с газом // Вычисл. технологии. 2002. Т. 7, № 3. С. 72-84.
[4] БисикАло Д.В., Боярчук А.А., КАйгородов П.В., Кузнецов О.А. Численное моделирование переноса вещества в тесных двойных звездах на компьютерах с параллельной архитектурой // Математическое моделирование: Проблемы и результаты. М.: Наука, 2003. С. 71-94.
[5] Попов М.В., Устюгов С.Д., Чечеткин В.М. Численное моделирование крупномасштабной неустойчивости при взрыве сверхновой II типа // Там же. С. 95-122.
[6] Четверушкин Б.Н., Тишкин В.Ф. Применение высокопроизводительных многопроцессорных вычислений в газовой динамике // Там же. С. 123-168.
[7] АнтонЕнко М.Н., Конюхов А.В., КрАгинский Л.М., Опарин А.М., Форто-вА С.В. Применение универсальной технологии параллельных вычислений в пространственных задачах физической газодинамики и сейсмического моделирования // Там же. С. 169-198.
[8] Четверушкин Б.Н. Кинетически согласованные схемы в газовой динамике. М.: Изд-во МГУ, 1999.
[9] Шильников Е.В., Шумков М.А. Моделирование трехмерных нестационарных течений газа на МВС с распределенной памятью // Мат. моделирование. 2001. Т. 13, № 4. C. 35-46.
[10] Balls G.T., Baden S.B., Colella P. SCALLOP: A Highly Scalable Parallel Poisson Solver in Three Dimensions // Supercomputing, 2003.
[11] Norman M.L. Introducing ZEUS-MP: A 3D, parallel, multiphysics code for astrophysical fluid dynamics // Astrophysical Plasmas: Codes, Models and Observations, eds J. Arthur, N. Brickhouse and J. Franco, Rev. Mex. Astron. Astrophys. (Conf. Ser.). 2000. Vol. 9. P. 66-71.
[12] Springel V. The cosmological simulation code GADGET-2 // Mon. Not. Roy. Astron. Soc. 2005. Vol. 364, N 4. P. 1105-1134.
[13] Nipoti C. Dissipationless collapse, weak homology and central cores of elliptical galaxies // http://arxiv.org/astro-ph/0605260
[14] Shirokoy A., Bertschinger E. GRACOS: Scalable and load balanced P3M cosmological N-body code // http://arxiv.org/astro-ph/0505087
[15] Genovese L., Deutsch T., Neelov A., Goedecker S., Beylkin G. Efficient solution of Poisson's equation with free boundary conditions // http://arxiv.org/cond-mat/0605371
[16] Четверушкин Б.Н. Проблемы эффективного использования многопроцессорных вычислительных систем // Информац. технологии и вычисл. системы. 2000. № 2. C. 22-34.
[17] Kochevsky A.N. Possibilities for simulation of fluid flows using the modern CFD software tools // http://arxiv.org/abs/physics/0409104
[18] Frenk C.S. et AL. The Santa Barbara Cluster comparison project: A comparison of cosmological hydrodynamics solutions // Astrophys. J. 1999. Vol. 525, N 2. P. 554-582.
[19] ПолячЕнко В.Л., Фридман А.М. Равновесие и устойчивость гравитирующих систем, М.: Наука, 1976.
[20] Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике. М: Наука, 1982.
[21] Toro E.F. A linearised Riemann Solver for the time dependent Euler equations of gas dynamics. Proc. Ray Soc. L. 1991. Vol. A434. P. 683-693.
[22] Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Мат. сб. 1959. Т. 47, вып. 3.
[23] Куликов И.М. Численное моделирование вращения газа в гравитационном поле // Тр. XLIII Междунар. науч. студ. конф. "Студент и научно-технический прогресс". Новосибирск: НГУ, 2005. C. 207-212.
[24] Вшивков В.А., Лазарева Г.Г., Куликов И.М. Операторный подход для численного моделирования гравитационных задач газовой динамики // Вычисл. технологии. 2006. Т. 11, № 3. С. 27-35.
[25] Яненко Н.Н. Метод дробных шагов многомерных задач математической физики. Новосибирск: Наука, СО, 1967.
[26] Kraeva M.A., Malyshkin V.E. Assembly technology for parallel realization of numerical models on MIMD-multicomputers // Future Generation Comp. Syst. 2001. Vol. 17, N 6. P. 755-765.
[27] Вшивков В.А., Тарнавский Г.А., Неупокоев Е.В. Распараллеливание алгоритмов прогонки: многоцелевые вычислительные эксперименты // Автометрия. 2002. Т. 38, № 4. С. 74-86.
[28] http://www.fftw.org
Поступила в редакцию 11 сентября 2006 г., в переработанном виде — 1 февраля 2007 г.
