CC BY
47
УДК 123.456, 789.012
ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ РЕШЕНИЯ ЗАДАЧ ДИНАМИКИ ЗАРЯЖЕННЫХ ЧАСТИЦ С УЧЕТОМ БАЛАНСИРОВКИ ВЫЧИСЛИТЕЛЬНОЙ НАГРУЗКИ1
Е.А. Берендеев, М.А. Боронина, В.Д. Корнеев
Рассмотрены задачи динамики встречных пучков заряженных частиц в ускорителях и динамики плазменных электронов в ловушке с инверсными магнитными пробками и муль-типольными магнитными стенками. Модели построены на основе метода частиц в ячейках. Такие задачи требуют большого объема вычислений и могут быть решены только с применением мощных суперЭВМ. Для равномерной и полной загрузки вычислительных узлов выполнена модификация эйлерово-лагранжевой декомпозиции в случае существенно неравномерного распределения частиц по пространству и по времени.
Ключевые слова: метод частиц в ячейках, параллельные алгоритмы, балансировка нагрузки, физика плазмы, встречные пучки, ускорители частиц.
Введение
Данная работа посвящена созданию и исследованию параллельных алгоритмов для моделирования некоторых задач физики плазмы методом частиц.
В частности, для моделирования динамики встречных пучков заряженных частиц в современных ускорителях с учетом высоких значений релятивистского фактора частиц и трехмерности задачи. При этом требуется проводить расчеты на сетках порядка 512x512x512 с количеством частиц в пучке ~1010. В основе существующих стандартных численных моделей и алгоритмов, направленных на решение задач с высокими релятивистскими факторами (103 и более), лежит разделение пучка вдоль оси коллективного движения на слои частиц. При этом взаимодействие встречных пучков сводится к попарному взаимодействию слоев: частицы одного слоя через поле поперечных сил влияют на динамику частиц другого. К настоящему времени созданы и параллельные алгоритмы решения задач динамики заряженных частиц в ускорителях. Обычно каждый процессор отвечает за свой слой частиц или за несколько слоев, в некоторых случаях имеется динамическая балансировка загрузки процессоров. Но используемый квази-трехмерный подход даже с учетом распараллеливания затрудняет учет продольных эффектов при критически высокой плотности пучка, когда за короткий промежуток времени пучок деформируется или даже разрушается, а также возникают трудности моделирования таким способом при сравнительно большом угле пересечения пучков (~20 мрад). Существенно нелинейное распределение плотности пучка (по гауссу с условием фокусировки), значительное изменение формы пучка при пролете через область (форма песочных часов), требования к количеству частиц пучка при заданной трехмерной сетке [1] приводят к необходимости создания кода с балансировкой нагрузки процессоров.
1 Статья рекомендована к публикации программным комитетом научной конференции Между-
народной научной конференции «Параллельные вычислительные технологии - 2014».
Аналогичные требования к размеру сетки и количеству частиц предъявляет задача динамики плазменных электронов в ловушке с инверсными магнитными пробками и мультипольными магнитными стенками. Величина магнитного поля внутри и на торцах ловушки существенно отличается. При этом, чтобы точно воспроизвести движение частиц в сильном магнитном поле, необходимо использовать достаточно мелкий временной шаг. В слабом магнитном поле величину временного шага можно увеличить и тем самым ускорить вычисления, т.к. траектории частиц в этом случае достаточно гладкие. Однако при построении параллельного алгоритма необходимо учитывать, что за реальный промежуток времени разные частицы совершают движение за разное количество временных шагов. Таким образом, объем вычислений для каждой частицы будет разный, даже если частиц на каждом процессоре поровну.
Рассмотренные задачи обладают существенной неоднородностью распределения частиц как по пространству, так и по времени. Поэтому нашей целью явилось создание параллельного алгоритма с балансировкой нагрузки на процессоры, обладающего хорошей масштабируемостью для возможности проведения численных экспериментов с использованием высокого количества частиц. В работе описана постановка задачи, используемые модели и методы, а также описание специфики распараллеливания и сравнение с другими работами мирового уровня. Представлена блок-схема алгоритма и результаты исследований его для обеих задач: время счета, эффективность распараллеливания, ускорение. Приведены примеры использования алгоритма для проведения математического моделирования соответствующих физических явлений.
1. Постановка задачи и методы решения
В задаче рассматривается движение частиц, которое происходит в вакууме в самосогласованных электромагнитных полях с учетом внешней составляющей электромагнитного поля. При этом модельный пучок представляет собой набор достаточно большого количества частиц с соответствующими параметрами [2, 3]. Такое движение моделируется с помощью кинетического уравнения Власова [4] для функции распределения
частиц с положительным f+ = f+ (г, р, t) и отрицательным зарядом /- = /- ( Г,Р,* ) :
д+- - д+- * дт
+у+- -
ді ’ дг
+ F+-^- = 0. др
(1)
Сила Лоренца , действующая на заряженную частицу, определяется из соотношения
,- = я+,-
\
(2)
Импульс частицы р + - связан со скоростью релятивистским фактором у+ _ равенством р+ _ = у+ _ ту+ _, при этом у+- = 1^1 — ^ С1 , где с - скорость света.
Система уравнений Максвелла, связывает между собой плотности заряда п+ , п , ток j и напряженности электрического и магнитного полей Е и Н :
ше=-1 д£,
с Э/
- 4п - 1 дЕ (3)
гоШ = —; +------------, (3)
с с д/
п ~+
divE = 4п (n_e + n_e+ ) divH = 0.
Входящие в эти уравнения плотность заряда и плотность тока определяются через моменты функции распределения частиц:
n+
=J /+dv.
V
= J f-dV,
V
J = J (f+v+ e+ + fy-Є)dV.
V
V
(4)
n
Уравнения характеристик уравнения Власова совпадают с уравнениями движения частиц:
-
dt (5)
іїг+_ -
~аТ~У+> - ‘
Система уравнений Власова-Максвелла решается методом частиц в ячейках, с использованием ядра PIC. Так как в уравнения входят первые производные и по времени, и по пространству, то применяется схема с перешагиванием [5].
Компоненты магнитного поля H вычисляются в центрах ребер ячеек, образованных пространственной сеткой, а значения электрического поля E вычисляются в серединах граней этих ячеек (рис. 1).
Рис. 1. Выбор узлов сетки для аппроксимации
При этом электрическое поле вычисляется в целые моменты времени пт, а магнитное поле - в дробные моменты времени п(т + 1/2). Токи вычисляются в тех же точках пространства, что и электрическое поле, но по времени берутся в дробных узлах. Значения импульсов частиц также вычисляются в моменты времени п(т + 1/2), а координаты
- в целые. При этом все производные, участвующие в уравнениях, как по времени, так и по пространству, записываются через центральные разности, что обеспечивает при использовании этой схемы второй порядок по времени и по пространству.
Сила Лоренца определяется линейной интерполяцией по каждому направлению сил электромагнитных полей, вычисленных в каждом из ближайших к частице восьми узлов.
Уравнения Максвелла решаются по следующим конечно-разностным схемам:
1 1
— т+— — т—
Н 2 - Н 2 -
___________Н____________-—гО/ Ет
Т ОТь (6)
77т+1 77 т 1 1
Е — Е — т+— — т+—
- у 2 + гоік Н 2
где
HZг, к, / - к - 1,/ - Ну і, к - , / — Ну і, к - , / - 1
к
к
НХі - к, / — НХі - к, / - 1 Н2і, к, / - ^ — Н2і - 1, к, / - ^
к
к
Ну і, к - 2, / — Ну і - 1, к - 2, / НХі - -2, к, / — НХі - у
к
к
77 1 1 17 1 1 17^ 1^1 1
Е2г -k + у,1 - Е2г -k -1 Еуг -k, I + - - Еуг -k, I - -
11 2’k - 2
2 ’ ’ 2
2 ’ ’ 2
К
к
Шн Е =
1 1 1 1 1 1 1 1
ЕХг, k - 2,1 + 2 — ЕХг, k - 2’ I - 2 Е21 + у, k - у, 1 — Е^г - у, k - у, 1
(7)
к
к
77 1 1^1 1 г 1 1 г 11
Еуг + 2, k, 1 - 2 — Еуг - у, k, 1 - У ЕХг, k + у, 1 - у — ЕХг, k - у, 1 - у
22
2 ’ 2
11
— I-
2 ,1 2
К
к
Импульсы каждой частицы находятся с помощью выражений:
= 4
Е т +
1 1
тл— т—
’ 2 + V 2
Нт
(8)
т
2
2
2
Токи вычисляются по координатам частиц с использованием ядра Р1С-метода по схеме, предложенной Вилласенором и Бунеманом [6]. Такой метод вычисления токов позволяет автоматически удовлетворить разностному уравнению неразрывности и, следовательно, точно выполнить разностный закон Гаусса. Это значительно уменьшает ошибки аппроксимации и делает алгоритм более устойчивым.
2. Параллельная реализация алгоритма
Метод частиц, созданный еще в 60-х годах 20 века, в настоящее время является развитым и широко применяемым методом для решения задач бесстолкновительной динамики частиц, ежегодные публикации демонстрируют постоянный интерес к этому методу.
В качестве приложения описанного выше алгоритма используются две задачи. Одна из них — получение мощных нейтральных пучков для установок управляемого термоядерного синтеза. Наиболее эффективным методом получения таких пучков является нейтрализация пучков отрицательных ионов в плазменной ловушке — мишени. В ИЯФ СО РАН предложена линейная ловушка с обратным магнитным полем [7]. Для ограничения радиальных потерь плазмы используются мультипольные магнитные стенки кольцевой геометрии. В осесимметричной ловушке с кольцевым магнитным полем отсутствует азимутальный компонент поля, а так же отсутствует стационарное азимутальное электрическое поле. Оценка и минимизация потерь плазмы в широко апертурные проходные отверстия в торцах, в которых находятся инверсные магнитные пробки, а также через цилиндрические мультипольные магнитные стенки ловушки на ее вакуумную камеру, может быть исследована только с помощью математического моделирования. При этом, наиболее полно динамика плазменных электронов может быть описана уравнением Больцмана
dt
dr
op
(9)
Здесь St{ f+-} — функция, описывающая следующие физические процессы:
• ионизация атома водорода,
• ионизация и диссоциация молекулы Н2
• диссоциативное возбуждение и диссоциативная рекомбинация Н2+
• диссоциативная рекомбинация D2+
• перезарядка протонов на атомах водорода.
Решение уравнения Больцмана (9) можно свести к решению уравнения Власова (1) и корректировке траекторий частиц с учетом рассеяния, используя методы Монте-
В настоящей работе рассеяние не учитывается, оно будет являться предметом дальнейшего исследования. Таким образом, решается система, описанная в параграфе 1, но в цилиндрической системе координат. Задача рассматривается в двумерной R-Z геометрии. Особенностью данной задачи является движение частиц под воздействием магнитных полей с большими градиентами. Для ускорения расчета траекторий частиц использовался динамический шаг по времени, обеспечивающий изменение напряженности магнитного поля не более 20% за один шаг.
Другой вариант использования алгоритма — трехмерное численное моделирование эффектов встречи пучков частиц с критическими параметрами (в частности, с высокими значениями релятивистских факторов и с предельно высокой плотностью зарядов) в ускорителях. В этом случае пучок может не только сильно сжиматься, но и разрушаться, поэтому требуется исследование устойчивости пучков. С точки зрения математического моделирования в задачах с большими значениями релятивистского фактора (у ~ 103) имеется существенное отличие от задач с малыми значениями (у ~ 5). Известно, что поле движущейся заряженной частицы в лабораторной системе координат вытягивается в у раз поперек оси движения и сокращается в уг раз вдоль этой оси. Так, например, при значениях релятивистского фактора у ~ 103 отношение поперечных размеров к продольному, на которых поля близки по абсолютной величине, составляет ~ 109, и использование традиционных путей решения становится невозможным.
В основе существующих стандартных численных моделей и алгоритмов, направленных на решение задач с высокими релятивистскими факторами (103 и более), лежит разделение пучка вдоль оси коллективного движения на слои частиц. При этом взаимодействие встречных пучков сводится к попарному взаимодействию слоев: частицы одного слоя через поле поперечных сил влияют на динамику частиц другого (коды Guinea-Pig, ODYSSEUS) [9, 10]. В большинстве существующих программных кодов, основанных на данном подходе, для получения поперечного поля используются либо формулы Бас-сетти - Ерскине [11], либо двумерное уравнение Пуассона с использованием быстрого преобразования Фурье (при этом граничные условия вычисляются через двумерную функцию Грина). К настоящему времени созданы и параллельные алгоритмы решения задач динамики заряженных частиц в ускорителях (коды COMBI, IMPACT, Beam Beam 3D) [12-13]. Обычно каждый процессор отвечает за свой слой частиц или за несколько слоев, в некоторых случаях имеется динамическая балансировка загрузки про-
цессоров. Но используемый квазитрехмерный подход затрудняет учет продольных эффектов при критически высокой плотности пучка, когда за короткий промежуток времени пучок деформируется или даже разрушается, а также возникают трудности моделирования таким способом при сравнительно большом угле пересечения пучков.
При этом интерес представляет изучение эффектов на расстоянии порядка дисперсии в поперечном направлении. Сама эта величина является небольшой, за счет фокусировки пучка (hour-glass эффект) его размеры в поперечном направлении сильно увеличиваются. Поэтому необходимо проводить исследования эффектов при больших размерах области на малых расстояниях от ее центра, то есть брать достаточно много и мелких пространственных шагов в поперечном направлении, что также приводит к необходимости использования параллельного алгоритма. Кроме того, временной шаг требуется выбирать не из соображений сходимости, так как он уже достаточно мал, а из соображений устойчивости метода. Таким образом, уменьшение сетки в поперечном направлении ведет не только к увеличению количества действий, связанных с количеством узлов сетки, но и к увеличению количества шагов программы.
Повышение точности расчетов методом части может быть достигнуто увеличением количества узлов пространственной сетки, однако такой способ не всегда возможно применять ввиду ограниченности памяти ЭВМ. Кроме того, для решений методом частиц характерны численные шумы, «самосила» и другие негативные эффекты, связанные с введением пространственной сетки. Одним из способов уменьшения влияния этих эффектов является увеличение количества частиц в ячейке. Кроме того, при увеличении количества узлов по каждому направлению, количество частиц также должно быть увеличено. Поэтому единственным вариантом достижения требуемых количественных характеристик является распараллеливание.
Поэтому для распараллеливания алгоритма для обеих задач выбран метод декомпозиции. Наиболее простая декомпозиция представляет собой разделение области на полосы в соответствующей системе координат. Для наибольшей эффективности вычислений с целью уменьшить число межпроцессорных коммуникаций необходимо модифицировать имеющиеся алгоритмы с учетом специфики решаемых задач [14]. Модификация представляет собой эйлеро-лагранжеву декомпозицию [15].
Рис.2. Эйлерово-лагранжева декомпозиция.
С каждой подобластью-полосой связана группа процессоров, частицы в каждой подобласти разделены между всеми процессорами группы, поэтому, чем больше частиц в определенной области - тем больше процессоров в соответствующей группе. Выбором количества процессоров в группе можно добиться необходимого максимального количества частиц в процессоре. Схема распараллеливания приведена на рисунке
2. Таким образом, в случае неравномерного распределения частиц по пространству, на каждую полосу приходится количество процессоров пропорциональное количеству частиц в полосе.
Если же частицы движутся с различным шагом по времени, то частицы распределяются между процессорами полосы не равномерно. Поскольку величина магнитного поля в пределах ячейки существенно не меняется, можно для каждой ячейки ввести поправочный коэффициент, характеризующий величину шага по времени для частиц данной ячейки. Чем больше магнитное поле - тем меньше поправочный коэффициент. Исходя из этих рассуждений, для того, чтобы вычислительная нагрузка на процессоры была сбалансирована, число процессоров Nт в группе оценивается следующим образом:
Здесь Nj — число частиц в ячейке j, tj — поправочный коэффициент, N — общее число процессоров.
На рисунке 3 приведена блок-схема алгоритма. Здесь белым цветом обозначены блоки программы, выполняемые каждым процессором, серым — нулевым процессором,
и темно-серым ---- главными процессорами каждой группы. Индексом g обозначены
значения функция, вычисленные каждым процессором группы, которые далее используются главным процессором группы для получения соответствующей суммарной функции. Обозначение send описывает пересылку значений каждым процессором группы главному процессору группы, а обозначение recv — принятие значений главным процессором группы от всех процессоров его группы, обозначение bcast — отправку значений главного процессора всем процессорами группы и соответствующее получение их. В случае ловушки требуются дополнительные процедуры по добавлению/удалению частиц. Этого не требуется в случае встречных пучков, но взамен необходимо вычислять граничные условия на каждом шаге по времени. Поэтому переход к каждому из блоков показан пунктирной линией.
После вычисления начальных условий выполняется цикл по времени до достижения нужного момента времени. В этом цикле требуется решать уравнения Максвелла, и каждая группа решает их только в своей подобласти. Уравнение Максвелла для вычисления электрического поля решается только главными процессорами группы, для этого производится сбор необходимых данных — плотностей и токов, а также обмен значениями на смежных гранях. За счет использования дополнительных узлов и равных значений электрического поля во всех процессорах магнитное поле вычисляется каждым процессором без пересылок. Также группы обмениваются частицами, оказавшимися в соответствующей подобласти. Такое распараллеливание по пространству и по частицам
(10)
позволяет решить проблему существенно неравномерного распределения плотности частиц в области, обеспечивая равномерную загрузку процессоров внутри группы.
- вычисление пробных координат - определение N3, в каждой группе - вычисление координат и других начальных значений 1
Получение начальных значений от 0-вого процессора
1 1 1
- вычисление pg0 - send pg0 - вычисление pg0 - send pg0 - вычисление pg0 - send pg0
1 1 1
- гест pg0, Р°=Х pg0 - вычисление Е0 - гecv pg0, Р°=Х Pg0 - вычисление Е0 - гест р80, р0=Х Pg0 - вычисление Е0
1 1 1
Обмен смежными значениями Е0 между главными процессорами групп
1 1 1
bcast Е0 bcast Е0 bcast Е0
І І 1
Вычисление Н^1/2 Вычисление ^-1/2 Вычисление Н^1/2
- вычисление Hgm - вычисление vm+1/2 - вычисление хт+1/2 - вычисление Hgm - вычисление vm+1/2 - вычисление х^2 - вычисление Hgm - вычисление vm+1/2 - вычисление х^2
1 1 1
Обмен параметрами частиц между процессорами соседних групп
1 1 1
- вычисление jgm+1/2 - send jgm+1/2 - вычисление jgm+1/2 - send jgm+1/2 - вычисление jgm+1/2 - send jgm+1/2
1 1 1
- gatheг jgm+1/2, jm+1/2=^jgm+1/2 - вычисление Hm+1/2 - gatheг jgm+1/2, jm+1/2=^jgm+1/2 - вычисление Hm+1/2 - gatheг jgm+1/2, jm+1/2=£jgm+1/2 - вычисление Hm+1/2
добавление /удаление частиц (ионизация / рекомбинация) добавление /удаление частиц (ионизация / рекомбинация) добавление /удаление частиц (ионизация / рекомбинация)
л +
К Є
Е м 53 й
''5 &
а ^ >к к я о Ч о
0>
К
К
й
£
СП
Конец расчета
Рис. 3. Блок-схема алгоритма
3. Эффективность параллельного алгоритма
Описанный алгоритм был протестирован на некоторых характерных примерах движения ультрарелятивистского пучка. В частности, предполагается, что пучок моноэнер-гетических электронов двигается вдоль оси z. Плотность частиц распределена по закону Гаусса со следующими параметрами фокусировки в безразмерных величинах Sx = Sy =
_7
540" , Px = Py=0,1, Gz=0,1, где Sx and Sy горизонтальный и вертикальный эмиттансы
пучка, Px, Py — соответствующие значения бета-функции, Gz — размер пучка по оси z.
Релятивистский фактор у = 6,85-103 , заряд пучка Q= 2,63408e. Размер области в без-
_2
размерных величинах Lx = Ly = 10 , Lz = 1. Размер сетки 120х120х120 узлов, временной шаг 10-5, количество шагов 3-104.
Все численные эксперименты проводились на кластере Сибирского Суперкомпью-терного Центра (ИВМиМГ), с 576-ю 4-ядерными процессорами Intel Xeon E5450/E5540/X5670.
На рис. 4 показана зависимость поля Ex от координаты z в плоскости (x,z) при раз-
5 б 7
ном количестве частиц 10 , 10 , 10 в конечный момент времени в безразмерных величинах.
Рис. 4. Зависимость поля Ех от координаты z в плоскости (х^) во всей области (слева) и в более мелком масштабе (справа).
Из рисунков видно, что при увеличении количества частиц остаточное поле, следующее за пучком, уменьшается, сглаживается. Рисунки также демонстрируют сходимость метода. Если же количество частиц в ячейке недостаточно, то высокие значения шумов в методе частиц приводят к неустойчивостям, затрудняющим исследование физических эффектов.
В табл. 1 представлены время расчета пролета пучка через область и максимальное количество частиц в процессоре, как в начальный момент (предпоследняя колонка), так и в течение всего расчета (последняя колонка). Таблица демонстрирует преимущество добавочных процессоров в группах, соответствующих высокой плотности частиц. Для линейной декомпозиции с количеством процессоров Np=Npg=12 время вычислений составило 84 с, и 98 с — для Np=Npg=6, но время расчета намного меньше, когда используется всего 12 процессоров в 6 группах. Таким образом, намного эффективнее
Таблица 1
Время расчетов и максимум частиц в процессоре
N«5 500т, с jmax, 106, 1 = 0 • 1 |-|6 jmax, 10
6 6 98 475025 500162
6 10 39 159123 167167
10 10 93 396014 497950
10 12 49 198030 248795
12 12 84 120590 494225
6 12 37 119283 125476
добавить 4-6 процессоров в центральные группы и получить время 37-39 с, чем распределить эти процессоры в линию и получить время 93 с. Этот эффект связан с существенным уменьшением числа частиц (с 4-105 до 1-1,5-105) в начальный момент времени, что влечет за собой и уменьшение количества частиц по всему расчету до 1,2—1,8-10 . Различие между 10 и 12 процессорами в 6 группах невелико вследствие межпроцессорных пересылок внутри группы.
Рисунок 5 демонстрирует зависимость эффективности распаралеливания от количества процессоров N (^) для Npg=5, сетки 60х60х60 и числа частиц 106. Время счета 1000 временных шагов составило 300 с. Эффективность распараллеливания рассчитывалась по формуле Eff = ЩТы/Тт, где N — количество процессоров, N0 - начальное количество процессоров. В идеальном случае необходимо сравнивать с монопроцессорным вариантом алгоритма, N0=1, но для проведения такого расчета недостаточно ресурсов памяти ЭВМ. Поэтому мы считали, что эффективность равна 1,0, когда мы используем N=Npg=5. Тм и Тш — время расчетов с соответствующим количеством процессоров.
Рис. 5. Эффективность параллелизации Рис. 6. Эффективность параллелизации
^=5, J=106, сетка 60x60x60 N =20, J=107, сетка 120x120x60
Добавление новых процессоров в области высокой плотности существенно увеличивает эффективность метода, когда число процессоров не очень велико (~20 в 5 группах),
но далее обмен трехмерными массивами сеточных значений занимает все больше времени, а поэтому эффективность распараллеливания падает.
Аналогичный результат получен и для параметров N =20 на сетке 120х120х60, число частиц 107, время расчета 1000 шагов по времени 67 с (рис. 6.). При этом провалы на графике связаны с несимметричным распределением процессоров в группах, поэтому некоторых из них простаивают. Из численных экспериментов следует, что оптимальное соотношение процессоров Npg и N достигается, когда число групп Npg примерно в 7 раз меньше размера сетки вдоль направления распараллеливания, а общее число процессоров N больше количества групп Npg примерно в 3-5 раз. В случае большой сетки эффективность падает за счет пересылок копий сеточных значений, тем не менее - это единственный способ проведения расчетов с таким большим количеством частиц и с таким мелким шагом по времени. Т.к. разумный минимум узлов в процессоре N равен 4 (2 узла + 2 вспомогательных узла), то параметры расчета ограничены только ресурсами ЭВМ для массивов размером ~4^^, соответствующих сеточным значениям подобласти.
В качестве второй задачи рассматривается динамика плазменных электронов в ловушке мишени. Поскольку задача двумерна, объем пересылаемых данных существенно ниже. В этом случае возможно добиться достаточно высокой масштабируемости алгоритма.
Характеристики задачи: температура плазмы 5 эВ, размер области 6,1 см х 1,2 см. Сетка 4096x128 узлов, общее число модельных частиц 5 242 880 000. Расчеты проводились на суперкомпьютере «Ломоносов» с использованием до 8192 процессорных ядер.
В таблице 2 представлено время расчета одного шага (в секундах) при использовании различного количества процессорных ядер, а также полученное при этом ускорение. В связи с большим объемом требуемой оперативной памяти, масштабируемость рассматривается относительно 1024 процессорных ядер.
Таблица 2
Время счета одного шага для различного числа процессоров и полученное ускорение
Общее количество используемых процессорных ядер Число процессорных ядер на подобласть Время расчета одного шага, с. Ускорение по сравнению с 1024 процессорными ядрами
1024 16 1,411 1
2048 32 0,778 1,81
4096 64 0,362 3,87
8192 128 0,183 7,71
На рис. 7 показаны траектории движения плазменных электронов, имеющих энергия 5 эВ, а также траектории движения электронов, стартующих с катода с энергией 150 эВ. Рисунок показывает, что инверсные магнитные пробки на торцах достаточно хорошо удерживают ловушки и границе инверсных пробок. Для количественного описания потерь плазму в ловушке, в то же время присутствуют
И (тт)
120
30
60
90
0
(
Ъ (тт)
Рис. 7. Траектории движения электронов мишенной плазмы под воздействием магнитного поля
потери плазмы на стенках плазмы в дальнейшем будет учтено влияние рассеяния электронов плазмы.
Заключение
Реализован параллельный алгоритм решения задач динамики плазмы, основанный на модификации эйлерово-лагранжевой декомпозиции в случае неравномерного распределения частиц по пространству и времени. В качестве примеров применения разработанного алгоритма были рассмотрены две задачи - динамики встречных ультрареляти-вистских пучков заряженных частиц в ускорителях и динамики плазменных электронов в ловушке с инверсными магнитными пробками и мультипольными магнитными стенками. Для задачи динамики встречных пучков уменьшение остаточного поля, следующего за пучком, с увеличением счетных параметров подтвердило сходимость метода. Для задачи ловушки получены траектории движения плазменных электронов и точки выхода плазмы на стенки ловушки. Применение алгоритма в обоих случаях позволило достичь высокой масштабируемости и увеличить число частиц в расчетах до 1010. Использование большого количества частиц даст возможность проводить расчеты с другими конфигурациями полей и частиц, и достигать большей точности решений.
Работа выполнена при поддержке грантов РФФИ 14-01-00392, 14-01-31220, а также МИП СО РАН №105 и №130.
Литература
1. Вшивков, В. А. Трехмерное моделирование динамики ультрарелятивистских пучков заряженных частиц: особенности вычисления начальных и граничных условий / В.А. Вшивков, М.А. Боронина // Математическое моделирование. — 2012. — Том 24, №2. — С. 67-83.
2. Хокни, Р. Численное моделирование методом частиц. / Хокни Р., Иствуд Дж. — М.: Мир, 1962
3. Березин, Ю.А. Метод частиц в динамике разреженной плазмы. / Березин Ю.А., Вшивков В.А. — Новосибирск, «Наука», 1980.
4. Власов, А.А. Теория многих частиц. / Власов А.А. — М.-Л., ГИТТЛ, 1950, 348 с.
5. Langdon, A.B. Electromagnetic and relativistic plasma simulation models / Langdon A.B, Lasinski B.F. / / Meth. Comput. Phys. — 1976. — Vol. 16. — P. 327-366.
6. Villasenor, J. Rigorous charge conservation for local electromagnetic field solver / Villasenor J., Buneman O. // Computer Phys. Comm. — 1992. — Vol. 69. — P. 306-316.
7. Dimov, G.I. Feasible scenario of startup and burnup of fusion plasma in ambipolar D-T reactor / Dimov G.I. // Transactions of Fusion Science and Technology. — 2011. — Vol. 59, No.1T. — P. 208-210.
8. Birdsall, C.K. Particle-in-Cell Charged-Particle Simulation Plus Monte Carlo Collisions With Neutral Atoms, PIC-MCC / Birdsall C.K. // IEEE Trans. Plasma Sci. — 1991. Vol. 19, No. 2. — P. 65-83.
9. C. Rimbault. GUINEA-PIG: A tool for beam-beam effect study // EUROTeV workshop.
— 2006. — Vol. Daresbury, 26-27 April.
10. Anderson, E. B. ODYSSEUS: A Dynamic Strong-Strong Beam-Beam Simulation for Storage Rings / E. B. Anderson, T. I. Banks, J. T. Rogers / / International Computational Accelerator Physics Conference. — 1998.
11. Bassetti, M. Closed Expression for the Electric Field of a Two-Dimensional Gaussian Charge / M. Bassetti, G. Erskine // CERNISR-ISR-TH/80-06. — 1980.
12. Kabel, A. A Multi-bunch, Three-dimensional, Strong-strong Beam-beam Simulation Code for Parallel Computers / A. Kabel, Y. Cai. / / 9th European Particle Accelerator Conference. — 2004.
13. Qiang, Ji. Parallel Strong-Strong/Strong-Weak Simulations of Beam-Beam Interaction in Hadron Accelerators / Ji Qiang, Miguel Furman, Robert D. Ryne, Wolfram Fischer, Ta-naji Sen, Meiqin Xiao / / AIP Conference Proceedings. — 2003. — Vol. 693. — P. 278281.
14. Андрианов, А.Н. Подход к параллельной реализации метода частиц в ячейках / Андрианов А.Н., Ефимкин К.Н. // Препринты ИПМ — Москва, 2009 г. — №009, 20c.
15. Берендеев, Е.А. Реализация эффективных параллельных вычислений при моделировании больших задач физики плазмы методом частиц в ячейках / Берендеев Е.А., Ефимова А.А. // Мат. междунар. конф. «Параллельные вычислительные технологии». Новосибирск, 2012. — C. 380-385.
Берендеев Евгений Андреевич, аспирант, ИВМиМГ СО РАН (Новосибирск, Российская Федерация), evgeny.berendeev@gmail.com
Боронина Марина Андреевна, к.ф.-м.н., младший научный сотрудник, Институт вычислительной математики и математической геофизики (Новосибирск, Россия), boronina@parbz.sscc.ru.
Корнеев Владимир Дмитриевич, к.т.н., доцент, старший научный сотрудник, Институт вычислительной математики и математической геофизики (Новосибирск, Россия), korneev@sscc.sscc.ru.
Поступила в редакцию 26 февраля 2014 г.
Bulletin of the South Ural State University Series “Computational Mathematics and Software Engineering”
2014, vol. 3, no. 1, pp. 97112
PARALLEL ALGORITHM FOR SOLUTION OF PROBLEMS OF CHARGED PARTICLE DYNAMICS BY THE USE OF LOAD BALANCE
E.A. Berendeev, Institute of Computational Mathematics and Mathematical Geophysics
SB RAS (Novosibirsk, Russian Federation)
M.A. Boronina, Institute of Computational Mathematics and Mathematical Geophysics
SB RAS (Novosibirsk, Russian Federation)
V.D. Korneev, Institute of Computational Mathematics and Mathematical Geophysics
SB RAS (Novosibirsk, Russian Federation)
We present new solution methods for dynamics of counter charged beams in accelerators and plasma electrons dynamics in a trap with inverse magnetic mirrors and multipole magnetic walls. The models are based on the particle-in-cell method. These problems require extremely large computations and can be solved by the use of powerful supercomputers only. A modification of the euler-lagrangian decomposition is implemented in order to achieve full and balanced load of computational nodes in case of highly nonuniform particle distribution in space and time.
Keywords: particle-in-cell method, parallel algorithms, load balance, plasma physics, counter beams, particle accelerators.
References
1. Vshivkov V.A., Boronina M.A. Trekhmernoe modelirovanie dinamiki ultrarelyativ-istskikh puchkov zaryazhennykh chastic: osobennosti vychisleniya nachalnykh i granich-nykh uslovij. Mathematical Models and Computer Simulations, 2012, vol. 24, no 2, pp. 67-83.
2. Hockney R.W., Eastwood J.W. Chislennoe modelirovanie metodom chastits. Moscow: Mir, 1962.
3. Berezin Ju.A., Vshivkov V.A. Metod chastits v dinamike razrezhennoj plazmy. Novosibirsk, «Nauka», 1980.
4. Vlasov А.А. Teoria mnogikh chastits — Moscow-Leningrad., GITTL, 1950, 348 p.
5. Langdon A.B, Lasinski B.F. Electromagnetic and relativistic plasma simulation models. Meth. Comput. Phys., 1976, vol. 16, pp. 327-366.
6. Villasenor J., Buneman O. Rigorous charge conservation for local electromagnetic field solver. Computer Phys. Comm., 1992, vol. 69, pp. 306-316.
7. Dimov G.I. Feasible scenario of startup and burnup of fusion plasma in ambipolar D-T reactor. Transactions of Fusion Science and Technology, 2011, vol. 59, no. 1T, pp. 208210.
8. Birdsall C.K. Particle-in-Cell Charged-Particle Simulation Plus Monte Carlo Collisions With Neutral Atoms, PIC-MCC. IEEE Trans. Plasma Sci., 1991, vol. 19, no. 2, pp. 65-83.
9. Rimbault. C. GUINEA-PIG: A tool for beam-beam effect study. EUROTeV workshop, 2006, vol. Daresbury, 26-27 April.
10. Anderson E.B., Banks T.I., Rogers J.T. ODYSSEUS: A Dynamic Strong-Strong Beam-Beam Simulation for Storage Rings. International Computational Accelerator Physics Conference, 1998.
11. Bassetti M., Erskine G. Closed Expression for the Electric Field of a Two-Dimensional Gaussian Charge. CERNISR-ISR-TH/80-06, 1980.
12. Kabel A., Cai Y. A Multi-bunch, Three-dimensional, Strong-strong Beam-beam Simulation Code for Parallel Computers, 9th European Particle Accelerator Conference, 2004.
13. Qiang J., Furman M., Ryne R.D., Fischer W., Sen T., Xiao M. Parallel Strong-Strong/Strong-Weak Simulations of Beam-Beam Interaction in Hadron Accelerators. AIP Conference Proceedings, 2003, vol. 693, pp. 278-281.
14. Andrianov A.N., Efimkin K.N. Podkhod k parallelnoj realizatsii metoda chastits v yachejkakh. Preprinty IPM, Moscow, 2009, no. 009, 20 p.
15. Berendeev E.A., Efimova A.A. Realizatsija effektivnykh parallelnykh vychislenij pri modelirovanii bolshikh zadach fiziki plazmy metodom chastits v jachejkakh. Materialy mezhdunarodnoy konferentsii “Parallelnye vychislitelnye tekhnologii”. Novosibirsk, 2012, pp. 380-385.
Received 26 February 2014
