Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 3 (2), с. 285-292
УДК 519.63
СРАВНИТЕЛЬНЫЙ АНАЛИЗ ПРОИЗВОДИТЕЛЬНОСТИ КЛАСТЕРНЫХ СУПЕРЭВМ НА ПРИМЕРЕ ЗАДАЧИ О РЕЛАКСАЦИИ ЭЛЕКТРОННОГО ПУЧКА В ВЫСОКОТЕМПЕРАТУРНОЙ ПЛАЗМЕ1
© 2011 г. А.В. Снытников
Институт вычислительной математики и математической геофизики СО РАН, Новосибирск
Поступила в редакцию 20.01.2011
Проведено сравнение отдельных показателей производительности четырех отечественных кластеров, достигнутых на задаче о моделировании взаимодействия электронного пучка с плазмой. Модель построена на основе метода частиц в ячейках. По результатам сравнения можно сделать вывод, что важнейшим показателем быстродействия кластера на реальных задачах является скорость доступа к оперативной памяти.
Ключевые слова: параллельное программирование, метод частиц в ячейках, плазма, теплопроводность.
Актуальность работы
Актуальность настоящей работы связана с необходимостью повышения эффективности использования кластерных суперЭВМ. Для этого требуется определить показатели их быстродействия при решении реальных физических задач. Причина того, что нельзя ограничиться измерениями только на универсальных тестах типа ЫпРаек заключается в том, что существует очень большая разница между декларируемой производительностью той части кластера, которая доступна конкретному пользователю в конкретном расчете, и реально достигнутой. Например, если используется одна четвертая часть процессоров кластера с пиковой производительностью 5.4 терафлопс, то не будет большой ошибкой сказать, что пиковая производительность используемой части кластера должна быть на уровне 1 терафлопс. Фактически на рассматриваемой задаче получается 0.18 терафлопс, что вызвано не только недостатками программы пользователя, но и неприспособленностью кластера к такому типу задач.
Также возможно использование результатов данной работы при определении рейтинга кластеров. В настоящий момент списки Тор50 и Тор500 выстроены в порядке убывания пиковой
1 Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2010».
производительности и производительности на тесте ЫпРаек что, разумеется, дает определенную информацию о сравнительной скорости работы представленных там машин. Но очень многие факторы, такие как скорость работы и объем дисков, пропускная способность шины памяти и коммуникационной сети, неоднородность оборудования и т.д., - остаются за пределами рассмотрения. А это именно те проблемы, с которыми придется столкнуться при попытке посчитать на кластере большую задачу.
Физическая задача, на примере решения которой будет рассматриваться производительность кластеров, представляет собой следующее. На многопробочной магнитной ловушке ГОЛ-3 (ИЯФ СО РАН) наблюдается понижение электронной теплопроводности на 2-3 порядка по сравнению с классическим значением [1] в результате релаксации в плазме мощного электронного пучка. Этот эффект представляет большой интерес с точки зрения построения энергетического термоядерного реактора, и на данный момент теоретическое описание его отсутствует.
Необходимость применения суперкомпью-терных вычислений обусловлена тем, что требуется, во-первых, иметь достаточно подробную сетку для того, чтобы воспроизвести резонансное взаимодействие релятивистского электронного пучка с плазмой, и, во-вторых, большое количество модельных частиц для того, чтобы промоделировать возникающую в дальнейшем турбулентность.
Описание модели
Постановка задачи
Численная модель, используемая для решения задачи о релаксации пучка, состоит из уравнений Власова для электронной и ионной компонент плазмы и системы уравнений Максвелла [1-5]. В общепринятых обозначениях эта система имеет следующий вид:
д!\,е , 5Ле , ^ дДе ° ^ ГЕ . 1 г В]
^Т+ ° ^ = %гI Е +-[у>В]
д? от др ^ с
^ 4п . 1 ЭБ
rot B = — 1 +--------------------,
c c dt '
rot Б = -
1 dB
c dt ^
(1)
div Б = 4пр, div B = 0.
В данной работе используется алгоритм решения системы уравнений (1), описанный в работе [2]. Далее все уравнения будут приводиться в безразмерном виде. Для обезразмери-вания используются следующие базовые величины:
• скорость света с = 3*1010 см/с,
1 а14 -3
• плотность плазмы n0 = 10 см ,
• плазменная электронная частота юе = 5.6*104 n1/2 с-1.
Уравнения Власова решаются методом частиц в ячейках (PIC). В рамках этого метода решаются уравнения движения модельных частиц, которые являются уравнениями характеристик для уравнения Власова. Величины с индексом i соответствуют ионам, с индексом e -электронам.
dt
= -k(e + [v,B]),
dt
= vi,e, Pi,e = YV
Y
Здесь коэффициент к для электронов равен
1, а для ионов - отношению масс электрона и иона.
Для нахождения электрических и магнитных полей используется схема, в которой поля определяются из разностных аналогов законов Фарадея и Ампера [2]:
gm+1/2 gm-1/2
= _4njm+1/2 + rot h Bm+1/2.
Рассмотрим следующую постановку задачи. В начальный момент в трехмерной области решения (размер области Ь), которая имеет форму прямоугольного параллелепипеда:
0 < х < Ь, 0 < у < Ь, 0 < 7 < Ь,
находится плазма, состоящая из электронов и ионов. Модельные частицы распределены по области равномерно. Задается плотность плазмы и температура электронов, температура ионов считается нулевой. Дополнительно в области присутствуют электроны пучка, которые также распределены по области равномерно (предполагается, что пучок уже вошел в расчетную область). Электроны пучка отличаются от электронов плазмы тем, что они имеют кинетическую энергию направленного движения е=1 МэВ, а их температура равна нулю. Модельные частицы, соответствующие электронам пучка, имеют меньшую массу, нежели модельные частицы, соответствующие электронам плазмы (отношение их масс равно отношению плотности плазмы и плотности пучка).
Итак, исходными параметрами задачи являются: плотность и температура электронов плазмы, отношение плотности электронов плазмы к плотности электронов пучка, энергия электронов пучка:
1 а17 -3
• плотность электронов плазмы п0 = 10 см . Плотность плазмы намеренно завышена по сравнению с реальной плазмой на установке ГОЛ-3 для того, чтобы эффекты релаксации пучка проявлялись в течение более короткого времени;
• температура электронов плазмы Т0 = 1 КэВ;
• отношение плотности электронов пучка к плотности электронов плазмы а = 10-3;
• энергия электронов пучка е = 1 МэВ.
Параллельная реализация
Распараллеливание выполнено методом декомпозиции расчетной области по направлению, перпендикулярному направлению движения электронного пучка. Используется смешанная эйлерово-лагранжева декомпозиция. Сетка, на которой решаются уравнения Максвелла, разделена на одинаковые подобласти по одной из координат. С каждой подобластью связана группа процессоров (в том случае, когда вычисления производятся на многоядерных процессорах, процессором для единообразия будет именоваться отдельное ядро). Далее, модельные частицы каждой из подобластей разделяются между процессорами связанной с этой подобластью группы равномерно вне зависимости от координаты.
т
т
Рис. 1. Декомпозиция области. Область решения разбита вдоль координаты Y на 4 подобласти, частицы каждой подобласти равномерно распределены между четырьмя процессорами независимо от координаты. Различные символы, обозначающие частицы: круг, квадрат, ромб, звезда, означают принадлежность частиц к разным процессорам
Каждый из процессоров группы решает уравнения Максвелла во всей подобласти. Далее решаются уравнения движения модельных частиц. После этого происходит суммирование значений тока по всей подобласти. Один из процессоров группы производит обмен граничными значениями тока и полей с соседними подобластями и затем рассылает полученные граничные значения всем процессорам своей группы. В случае, если и уравнения Максвелла для подобласти, и уравнения движения всех частиц подобласти частиц решаются на одном процессоре, вычисления с частицами занимают в 10-20 раз больше времени.
Эффективность распараллеливания
120%
100%
80%
60%
40%
20%
256
N
512
1024
Рис. 2. Эффективность распараллеливания, вычисляемая по формуле (2) для небольшого числа процессорных ядер. Расчеты проведены на МВС-100К, МСЦ РАН
здесь Т1 - время счета с использованием N1 процессоров, Т2 - время счета с использованием N2 процессоров, £ь £2 - характерные размеры задач, в данном случае число узлов сетки по координате X. При этом размер задачи увеличивается пропорционально числу процессоров, т.е. нагрузка на отдельный процессор не возрастает. Цель такого определения эффективности - понять, насколько увеличивается время счета при увеличении числа процессоров и неизменной нагрузке на один процессор. В идеале время должно остаться тем же самым (к = 100% в идеале).
За основу для сравнения при этом берется модель, требующая большого времени вычислений (сетка 256х 128х 128 узлов, 150 частиц в ячейке). В расчетах по определению эффективности увеличивается только размер сетки по X, все остальные параметры остаются неизменными. Вопрос о том, насколько быстрее можно посчитать на суперкомпьютере задачу, которую можно посчитать и на настольной машине, в данной работе не рассматривается: такие задачи не представляют интереса с физической точки зрения.
На рис. 2 представлена эффективность для сравнительно небольшого числа процессорных
Рис. 3. Эффективность распараллеливания, вычисляемая по формуле (2) для большого числа процессорных ядер. Расчеты проведены на МВС-100К, МСЦ РАН
Основная цель создания параллельной программы для моделирования аномальной теплопроводности в плазме - возможность счета на больших сетках и с большим количеством модельных частиц. Поэтому эффективность распараллеливания вычислялась следующим образом
Т N 5 к = -2 х—2 х^2 х100% Т N
(2)
ядер, на рис. 3 приведена эффективность для большего числа ядер (> 1000). Из этих рисунков можно сделать два вывода: во-первых, падение производительности при переходе с 1024 ядер на 2048 меньше, чем при переходе с 256 на 1024, во-вторых, созданная программа способна эффективно использовать более тысячи процессорных ядер.
Сравнение производительности кластеров
Таблица 1
Основные характеристики кластеров, на которых производились расчеты, по материалам Российской суперкомпьютерной конференции «Научный сервис в сети Интернет: Масштабируемость, параллельность, эффективность»,
сентябрь 2009, www.supercomputers.ru
№ по списку Top50 Название Узлы Сеть Пиковая производительность, терафлопс
1 МВС-100К 4xXeon E5450, 3 GHz, 8.192 GB RAM Infiniband 4x DDR, InfiniBand/ Gigabit Ethernet 95.04
2 СКИФ-МГУ 2xXeon E5472, 3 GHz, 8.192 GB RAM InfiniBand/ Gigabit Ethernet/ СКИФ-ServNet + IPMI 60
14 СКИФ-Cyberia 2xXeon E5150, 3 GHz, 4 GB RAM QLogic InfiniPath/ Gigabit Ethernet/ СКИФ-ServNet 12.002
20 Кластер НГУ 2xXeon E5355, 2.66 GHz, 8.192 GB RAM Infiniband 4x DDR, Gigabit Ethernet / Gigabit Ethernet 5.4
Каждый временной шаг работы программы состоит из следующих действий:
• Расчет электрического и магнитного поля;
• Расчет движения модельных частиц;
• Вычисление новых значений плотности тока и заряда.
Дополнительно на отдельных временных шагах (как правило, каждый сотый шаг) проводится выдача диагностической информации, важнейшей частью которой являются фурье-образы основных величин (плотность заряда и тока, модуль электрического и магнитного поля).
Размер сетки во всех расчетах этого раздела 512х64х64, 150 модельных частиц в ячейке, число использованных процессорных ядер -160. Расчетная область поделена на 32 части вдоль координаты У, далее частицы каждой подобласти дополнительно поделены между 8 ядрами. Время работы всех этих процедур было измерено с помощью профилировщика §рго£ В каждом случае приводится время работы одного вызова процедуры, то есть общее время, затраченное на эту процедуру, поделенное на число вызовов.
Расчет движения модельных частиц
и скорость работы оперативной памяти
Для вычисления новых значений координаты и импульса модельной частицы используются значения электрического и магнитного поля. Каждая компонента поля хранится в отдельном трехмерном массиве. Таким образом, на каждом временном шаге для каждой модельной частицы происходит обращение к шести трехмерным
массивам. Модельные частицы расположены внутри расчетной области случайным образом. Если даже модельные частицы расположены рядом в массиве, где хранятся их координаты, то сами значения координат будут близкими только вначале. В дальнейшем модельные частицы перемешиваются.
Это означает, что обращения к трехмерным массивам, содержащим электрическое и магнитное поля, являются неупорядоченными, и использование кэш-памяти в данном случае не позволяет сократить время счета. Таким образом, время расчета движения модельных частиц в основном определяется именно быстродействием оперативной памяти (пропускной способностью шины памяти). На рис. 4 показано время, затраченное на вычисление движения частиц на одном временном шаге. Вывод о том, что время обращения к оперативной памяти является определяющим при счете с частицами, подтверждается сравнением времен, полученных для МВС-100К и СКИФ-МГУ. Процессоры, установленные на этих кластерах, являются близкими по своим характеристикам, в то время как время счета отличается почти в два раза.
Из рис. 4 также видно, что существуют большие возможности для оптимизации, например, путем сортировки частиц по координате, что позволит упорядочить обращения к памяти и использовать кэш. Так как расчет движения модельных частиц занимает большую часть времени (от 92% на СКИФ СуЬепа до 64% на СКИФ-МГУ) и это именно та часть программы, которая лучше всего распараллеливается, оптимизация этой процедуры может при-
вести к ухудшению ускорения и параллельной эффективности, однако уменьшение общего времени счета представляется более важной задачей.
Рис. 4. Время работы процедуры интегрирования уравнений движения модельных частиц на различных кластерных системах. Расчеты проведены на СКИФ Cyberia (ТГУ, Томск), МВС-100К (МСЦ РАН, Москва), СКИФ-МГУ (МГУ, Москва) и на кластере НГУ (НГУ, Новосибирск)
быстрое преобразование Фурье. Все локальные переменные этой процедуры полностью помещаются в кэш, поэтому на примере этой процедуры можно судить о том, каким может быть время «быстрого» счета, то есть с использованием только кэша, без обращений к оперативной памяти.
Пересылка модельных частиц и скорость работы межпроцессорных коммуникаций
Скорость работы межпроцессорных коммуникаций была измерена с помощью процедуры пересылки частиц между процессорами. Эта процедура включает в себя определение модельных частиц, вылетевших за пределы подобласти, принадлежащей данному процессору, и находящихся в буферном пограничном слое, перемещение этих частиц в буфер для пересылки и собственно пересылку в зависимости от номера процессора: если четный номер, то вначале «левому» процессору, то есть тому, чья подобласть расположена левее по координатной оси, а затем «правому». Если номер процессора нечетный, то наоборот: вначале «правому», потом «левому». Вместе с отправкой частиц на другой процессор происходит прием частиц, перелетевших в подобласть, принадлежащую данному процессору. Дополнительно перед пересылкой собственно частиц соседние процессоры обмениваются количеством частиц, которые необходимо переслать.
Рис. 5. Время работы процедуры, реализующей одномерное преобразование Фурье на различных кластерных системах. Расчеты проведены на СКИФ Cyberia (ТГУ, Томск), МВС-100К (мСЦ РАН, Москва), СКИФ-МГУ (МГУ, Москва) и на кластере НГУ (НГУ, Новосибирск)
Одномерное преобразования Фурье и скорость счета
Для того, чтобы отделить время счета от времени обращения к оперативной памяти, было рассмотрено время работы процедуры, реализующей одномерное преобразование Фурье (рис. 5). Это процедура fftc из библиотеки NAG, она принимает на вход одномерный комплексный массив размером 512 или 64 и вычисляет
Рис. 6. Время работы процедуры пересылки модельных частиц на различных кластерных системах. Расчеты проведены на СКИФ Cyberia (ТГУ, Томск), МВС-100К (МСЦ РАН, Москва), СКИФ-МГУ (МГУ, Москва) и на кластере НГУ (НГУ, Новосибирск)
Таким образом, время, показанное рис. 6, затрачивается на просмотр списка частиц и четыре пересылки типа МР1_$епЖесу. Количество перелетающих частиц не может быть большим
из физических соображений (то есть если возникают большие потоки частиц с процессора на процессор, это означает, что расчет физически некорректный), поэтому размер буфера для пересылки жестко задан - 5% массива частиц (на практике обычно еще меньше). Это означает, что объем пересылаемых данных порядка 10 Мб в каждой пересылке.
На рис. 6 видно, что наименьшее время на пересылки тратится на СКИФ СуЬепа и СКИФ-МГУ. Возможно, причиной этого является технология ServNet, используемая на машинах семейства СКИФ. Наибольшее время на пересылки израсходовано на МВС-100К, что можно объяснить большими масштабами и разнородностью этого вычислительного комплекса.
Выдача физической диагностики
Основной причиной падения эффективности распараллеливания в данной задаче является выдача физической диагностики: фурье-образы распределений поля, заряда, тока, скоростей и импульсов электронов и ионов, временные зависимости температуры, энергии частиц и полей, аномальный и классический коэффициент теплопроводности. Большое количество диагностики, а также однопроцессорный характер ее выдачи (сбор по всем процессорам и затем запись на диск одним выделенным процессором) являются причиной того, что именно за счет диагностик происходит потеря производительности программы.
С другой стороны, их нельзя исключить из рассмотрения и анализировать скорость работы только «вычислительной» части программы: формальный результат по ускорению и эффективности, безусловно, был бы лучше, но программа написана для получения физического результата, а, значит, для диагностики. Кроме того, целью настоящей работы является оценка фактической производительности кластеров именно с точки зрения конечного пользователя, - формальные показатели их производительности хорошо известны.
Для примера была выбрана выдача температуры электронов и ионов. Эта диагностика выдается на каждом временном шаге. Для вычисления температуры необходимо вычислить дисперсию импульсов частиц отдельно на каждом процессоре, сложить частичные суммы, и затем 0-й процессор выдает вычисленные значения в файл. Таким образом, показанное на рис. 7 время выдачи температуры включает в себя расчет, пересылку и выдачу в файл. Так же, как и в
случае с пересылкой частиц, лидируют представители семейства СКИФ. Отсюда можно сделать вывод, что основное время при выдаче диагностики расходуется именно на пересылку, и оптимизацию диагностики нужно будет начинать именно с этого.
Рис. 7. Время работы процедуры, выполняющей расчет температуры электронов и ионов на различных кластерных системах. Расчеты проведены на СКИФ Cyberia (ТГУ, Томск), МВС-100К (МСЦ РАН, Москва), СКИФ-МГУ (МГУ, Москва) и на кластере НГУ (НГУ, Новосибирск)
Вопрос об оптимизации диагностических сообщений является очень актуальным, но его решение возможно только в том случае, когда будет точно известно, какую информацию и в каком объеме нужно выдавать, а такая ситуация равнозначна решению физической задачи. В данный момент диагностики различных видов убираются и добавляются после каждого вычислительного эксперимента.
Для того чтобы понять, насколько сильно выдача диагностики нуждается в оптимизации, приведем конкретные цифры: на СКИФ СуЬегіа расчет температуры занимает 2.39% всего времени работы программы, на МВС-100К 4.17%, на СКИФ-МГУ 3.94% и на кластере НГУ 4.52%. Это самая времяемкая процедура диагностики, все остальные потребляют на порядок меньше времени. Соотношение времени работы диагностики и времени счета следующее: расчет движения частиц занимает намного больше времени (от 64% до 92%), расчет поля - намного меньше (от 1.6% на МВС-100К до менее чем
0.1% на СКИФ СуЬегіа).
Аномальная теплопроводность в вычислительных экспериментах
Актуальность настоящей работы связана с тем, что в экспериментах на многопробочной
магнитной ловушке ГОЛ-3 (ИЯФ СО РАН) наблюдается понижение электронной теплопроводности на 2-3 порядка по сравнению с классическим значением. Известно множество работ по моделированию теплопроводности в термоядерных установках. Вычисление температуры в этих работах производится в основном с помощью гидродинамических уравнений. Таким образом функция распределения электронов по энергиям предполагается максвелловской, что может не соответствовать действительности.
Для вычисления неравновесных распределений применяют, в частности, бессеточные модификации метода частиц. Однако для метода частиц в ячейках не известны критерии, позволяющие понять, насколько правильным является полученное в расчетах распределение температуры, и, таким образом, корректность моде-
лирования теплопроводности также оказывается под вопросом. В результате проведения вычислительных экспериментов выяснилось, что движение электронов пучка, первоначально равномерное и однонаправленное, после релаксации пучка становится вихревым. Это приводит к тому, что поток энергии электронов
2
-у
также приобретает сложную
структуру с модуляциями по осям X и У, и, более того, возникают подобласти, где поток энергии близок к нулю (рис. 8). На рисунке линия уровня с наименьшим значением (0.01) соответствует 1% начальной величины потока энергии электронов. Из рисунка видно, что в значительной части расчетной области поток энергии меньше начального.
Рис. 8. Линии уровня потока энергии электронов в плоскости XY, z = ZM/2, момент времени t = 91.7 (в единицах плазменного периода). Величина потока нормирована на его значение в момент времени t = 0
Заключение
Описана параллельная реализация модели взаимодействия электронного пучка с плазмой, приведены показатели ускорения и эффективности. Измерено и проанализировано время работы отдельных частей программы на различных кластерах. В результате выяснилось, что важнейшим показателем быстродействия кластера с точки зрения программы, реализующей метод частиц в ячейках, является скорость доступа к оперативной памяти, а основным сдерживающим фактором для эффективности распараллеливания является выдача физических диагностик. В отношении самой программы показано, что существуют большие возможности для ускорения счета.
Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований, гранты № 08-01-615 и № 08-01-622, а также интеграционных
проектов СО РАН № 103, № 113 и № 26 и гранта Президента Российской Федерации для государственной поддержки молодых российских ученых № МК-3562.2009.9.
Список литературы
1. Астрелин В.Т., Бурдаков А.В., Поступаев В.В. Подавление теплопроводности и генерация ионнозвуковых волн при нагреве плазмы электронным пучком // Физика плазмы. 1998. Т. 24. № 5. С. 45462.
2. Вшивков В.А., Вшивков К.В., Дудникова Г.И. Алгоритмы решения задачи взаимодествия лазерного импульса с плазмой // Вычислительные технологии. 2001. Т. 6. № 2.
3. Кролл Н., Трайвелпис А. Основы физики плазмы. М: Мир, 1975.
4. Григорьев Ю.Н., Вшивков В.А. Численные методы «частицы-в-ячейках», Новосибирск: Наука.
2000.
5. Бедсел Ч., Лэнгдон Б. Физика плазмы и математическое моделирование. М: Мир, 1989.
A COMPARATIVE PERFORMANCE ANALYSIS OF CLUSTER SUPERCOMPUTERS BY THE EXAMPLE OF ELECTRON BEAM RELAXATION IN HIGH TEMPERATURE PLASMA
A.V. Snytnikov
We compare some performance characteristics achieved by four domestic cluster supercomputers when solving the problem of simulation of electron beam interaction with plasma. The model is based on the Particle-in-Cell method. The results of this comparison clearly show that the most important cluster performance indicator for real physical problems is the RAM access speed.
Keywords: parallel programming, Particle-in-Cell method, plasma, heat conductivity.