CC BY
76
Сер. 10. 2012. Вып. 3
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
УДК 519.6: 621.38
И. П. Юдин, Е. Е. Перепелкин
ИСПОЛЬЗОВАНИЕ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ ПРИ ИССЛЕДОВАНИИ ПРОПУСКНОЙ СПОСОБНОСТИ КАНАЛА ТРАНСПОРТИРОВКИ ПУЧКА ИОНОВ С УЧЕТОМ ПРОСТРАНСТВЕННОГО ЗАРЯДА
1. Введение. Моделирование динамики пучка заряженных частиц является актуальной задачей ускорительной физики. Особенно она усложняется при учете эффекта пространственного заряда. Существуют различные подходы по моделированию таких процессов. К наиболее распространенным относится метод крупных частиц [1]. Однако для корректного учета влияния эффекта пространственного заряда на динамику пучка необходимо использовать большое (свыше 1 000 000) число крупных частиц (или макрочастиц). Существенно проблема усложняется при необходимости многократного пересчета модели, например с целью оптимизации параметров установки. В таких случаях происходит существенное увеличение времени моделирования. Подход, который изложен в настоящей работе, основан на применении параллельных вычислений на графических процессорах (GPU) с использованием среды программирования CUDA (Compute Unified Device Architecture). Для определенного класса вычислительных задач GPU дает производительность выше, чем центральный процессор CPU, от десятков до сотен раз. Показано, что задача моделирования динамики частиц с учетом пространственного заряда хорошо подходит для вычисления на GPU. Как результат применения массивно параллельных вычислений на GPU в среде программирования CUDA в п. 3 получено ускорение около 90 раз на вычислительном модуле Tesla C1060. Такой результат существенно сократил время вычислений, тем самым облегчив процесс оптимизации.
2. Математическая постановка задачи. Опишем математический аппарат для моделирования динамики пучка с учетом пространственного заряда в инжекционном канале синхротрона. Физически пучок состоит из большого количества частиц, которое может достигать величины порядка 1015. Конечно, такое количество частиц (будем называть такие (реально существующие) частицы микрочастицами) невозможно использовать при реальном моделировании на компьютере. Заметим, что для хранения в памяти компьютера 3-мерных координат частиц в этом случае потребуется 1015 х 3х 4 Байта = 12 ТБ. Здесь считается, что координата имеет одинарную точность float
Юдин Иван Павлович — ведущий научный сотрудник Лаборатории физики высоких энергий Объединенного института ядерных исследований (г. Дубна Московской обл.). Количество опубликованных работ: 210. Научные направления: математическое моделирование, численные методы, методы оптимизации, распределенные вычисления, ускорительная физика. E-mail: yudin@jinr.ru.
Перепелкин Евгений Евгеньевич — старший научный сотрудник Лаборатории физики высоких энергий Объединенного института ядерных исследований (г. Дубна Московской обл.). Количество опубликованных работ: 60. Научные направления: математическое моделирование, численные методы, методы оптимизации, распределенные вычисления, ускорительная физика. E-mail: pevgeni@jinr.ru.
© И. П. Юдин, Е. Е. Перепелкин, 2012
и требует соответственно в памяти 4 Байта. Понятно, что 12 ТБ оперативной памяти в настоящее время не реальная величина. Заметим, что это всего лишь координаты частиц, а при вычислении необходимо применять для каждой частицы и вектор скорости. Таким образом, требуется размер памяти порядка 24 ТБ. Также стоит отметить, что при использовании двойной точности double это приведет к удвоению размера, т. е. к 48 ТБ. Для выхода из такой ситуации были разработаны различные модели, описываемые ниже.
2а. Модель макрочастиц (Метод крупных частиц). Одним из подходов в описании пучка явилось представление его в виде, приемлемом для компьютерного моделирования числа макрочастиц. Каждая макрочастица является, по сути, объединением некоего количества микрочастиц. Так, например, если в исходном физическом пучке число микрочастиц составляет 1015, а при моделировании такому пучку ставится в соответствие пучок из 106 макрочастиц, то из этого следует, что каждая макрочастица содержит 109 микрочастиц. И заряд, и масса такой макрочастицы равны 106mmicro и 109qmicro соответственно, где mmicro и qmicro - масса и заряд микрочастицы.
В результате пучок представляется как набор таких параметров:
(ri, Pi,t), i = 1,...,N,
в котором N - число макрочастиц, ri и pi - вектор координат и импульса i-той макрочастицы в момент времени t. Отметим, что ri и pi, в свою очередь, могут быть определены, как усредненные значения по всем микрочастицам, содержащимся в i-той макрочастице.
Уравнением движения таких макрочастиц является второй закон Ньютона, который в рассматриваемой ситуации принимает вид
= г = 1 ,...,N, (1)
где Fi - сила, действующая на i-тую макрочастицу. Такая величина представима в форме суперпозиции следующих сил:
Fi = qi (Eext (ri, t) + Es (ri, t) + [Vi, Bext (ri)]), (2)
здесь Eext (ri,t) и Bext (ri) - внешние электрические и магнитные поля, действующие на пучок. Распределение таких полей должно быть известным (или вычисляться в процессе компьютерного моделирования). Отметим, что внешнее электрическое поле Eext (ri,t) в отличие от внешнего магнитного поля Bext (ri), помимо пространственной координатной зависимости, еще имеет временную зависимость. Поле Es (ri,t) является электрическим полем, созданным системой зарядов макрочастиц. Такое поле не задано и требует своего определения в каждый момент времени. Связано это с тем, что частицы в каждый момент времени изменяют свое положение и соответственно меняется распределение электрического поля Es (ri, t). Вопрос о Es (ri,t) не тривиален, и для его решения используются различные подходы. В купе данные подходы объединены под общим названием задачи об учете пространственного заряда, которая рассматривается ниже. Отметим: в рассматриваемой задаче предполагалось, что собственное магнитное поле пучка пренебрежимо мало по сравнению с внешним магнитным полем Bext (ri).
Для решения уравнения движения (1) необходимо дополнить его начальными и «граничными» условиями. Одним из вариантов таких условий является знание начальных координат (1) r(0) и скоростей v(0) макрочастиц в пучке. Это предположение
дает стандартное условие Коши
г4=0 = г(0), ^ = У(0), г = 1,...,М. (3)
Другой вариант - знание так называемого распределения макрочастиц «по ВЧ фазе». Такой подход используется в случаях, когда пучок начинает постепенно выходить из некоторой щели источника или известно, что каждая г-тая макрочастица пучка пересекла некую поверхность в некоторый свой момент времени Ь^. Математически такое начальное условие выглядит следующим образом:
= г(0), ^и = г = 1,...,М. (4)
Получается, что для каждой г-той макрочастицы задано еще время ^, в момент кото-
(0) (0) т-г
рого для макрочастицы известны ее положение гг и вектор скорости .По сути, условие (4) отличается от предыдущего (3) тем, что в первом все координаты и скорости известны в один и тот же момент времени для всех макрочастиц пучка. Во втором случае координаты и скорости известны для каждой г-той макрочастицы пучка только в соответствующие ей моменты времени ^.
Теперь зададим «краевые условия». Кавычки для слова «краевые условия» здесь употребляются в силу неполной эквивалентности используемого в математике понятия «краевое условие». Здесь под «краевыми условиями» будем понимать ограничения по координатам на положение частиц. Другими словами, существует некоторая область V С В3, в которой ищутся координаты частиц г^, г = 1,...,М, в каждый момент времени Ь. Если в результате решения уравнений движения получается, что 3 ] : г- <// V, тогда ]-тая макрочастица «исключается» из рассмотрения. То есть число макрочастиц, участвующих в движении пучка, уменьшается. Такие макрочастицы называются потерянными на структурных элементах установки. Следует отметить, что этот факт приводит к изменению распределения электрического поля Е8 (г,Ь). Граница области V, по сути, задает геометрию установки, тем самым определяя область, в которой может двигаться пучок. Считается, что область V задана. Однако стоит отметить, что для задачи оптимизации апертуры магнитной системы ускорителя требуется найти «наилучшую» область V. В такой постановке необходимо рассматривать, в некотором смысле, обратную задачу, о которой будет сказано далее.
Итак, приходим к следующей постановке задачи:
т4 (7гУг) = Чг (Ееж4 + Ея (п,г) + [лги Вехг (г*)]) ,
Ъ = а = *
Г= г(0\ V= ^(0\ г¿|í=0 = г(0), ^=0 = V((0), 1 < г < м, г е V,
(5)
где - масса покоя г-той макрочастицы; с - скорость света в вакууме. Как указывалось ранее, постановку (5) необходимо дополнить процедурой нахождения поля Ея (г,Ь). Таким образом, нашу задачу следует решать с учетом проблем с пространственным зарядом.
2б. Задача с пространственным зарядом.. Рассмотрим постановку задачи о поиске собственных полей пучка. Запишем систему уравнений Максвелла
д та _ п _1_ — Л V"/
д1 я ' моео '
шУЕ8
го1Ея
где р - плотность заряда в пучке; Es - собственное электрическое поле пучка, т. е. электрическое поле, которое создает заряд пучка с плотностью р; Bs - собственное магнитное поле пучка, которое появляется вследствие движения пучка. В случае, когда все частицы покоятся, плотность тока пучка Js = 0, а поле Es стационарно, т. е. Bs = 0. Поэтому для корректного учета эффекта пространственного заряда необходимо в выражение для силы (2) добавить вклад от Es и Bs пучка, т. е.
F = qi (Eext (r,t) + Es (ri,t) + [Vi, Bext (ri) + Bs (ri, t)]). (7)
Собственные поля могут быть найдены из уравнений (6), если к ним добавить граничные условия. Для оценки влияния магнитного поля на поперечную фокусировку пучка рассмотрим модель цилиндрического, продольно однородного непрерывного пучка. Такой пучок движется вдоль оси OZ. В этом случае излучением электрического поля можно пренебречь. В результате уравнения (6) примут вид
rotBs = M0J, divDs = р.
Учтем, что
J = pv, Ds = —е0Уф, Bs = rot A,
где ф - скалярный потенциал собственного электрического поля пучка, а A - векторный потенциал собственного магнитного поля пучка. Скорость частиц в пучке можно представить так:
v = ¡3cez + 5v.
Здесь предполагаем, что Sv < f3cez, т. е. продольная скорость во много раз выше поперечной. В этом случае можно найти распределение собственных полей
в
А = —фех
c
с
(8)
£о
В результате сила, действующая на пучок от собственных полей, принимает вид
Fs = q (Es + [v, Bs]) « q (^-Уф + f3cez, rot т. е.
Как видно из формулы (9), эффект собственного магнитного поля состоит в поперечной самофокусировке пучка. Выше было отмечено, что в решаемых в этой работе задачах можно пренебречь влиянием собственного магнитного поля Bs. Такое утверждение делается на основании того факта, что энергии пучков, которые рассматриваются в этой работе, удовлетворяют условию y2 < 1. В результате задача сводится только к отысканию собственного электрического поля Es. Возможны два основных способа нахождения Es. Первый - это метод PP (частица на частицу), второй - метод PIC (частица-в-ячейке). Опишем их особенности.
2в. Метод частица на частицу (Particle to Particle, PP). Данный подход основан на представлении электрического поля Es в виде суперпозиции полей макрочастиц, содержащихся в пучке, т. е.
1 N
4п£0 j= \п — Tj \
При вычислении по формуле (10) надо обращать внимание на величину |т — Tj \. Если она будет очень маленькой, то можно получить сингулярность в расчете электрического поля. Для решения этой проблемы используют представление макрочастицы в виде облака с радиусом R. В случае, когда облако г-той макрочастицы пересекается с облаком '-той макрочастицы, поле, которое действует на г-тую макрочастицу, вычисляется следующим образом:
Е« U ») = дз (ri " ' lri " rjl < R-
Увеличение радиуса облака приводит к ослаблению эффекта пространственного заряда, поэтому выбор оптимального размера облака представляет отдельную тему для исследования. Стоит отметить, что сильный рост электрического поля вблизи макрочастицы приводит к необходимости использования достаточно мелкого шага интегрирования по времени уравнений движения макрочастицы. Вычислительная емкость формулы (10) - порядка O(N2), что приводит к значительным затратам по времени. Для уменьшения количества операций следует применить следующий подход. Идея состоит в точном учете вклада электрического поля только от ближайших соседей, от остальных макрочастиц вклад рассчитывается по некой аппроксимационной формуле.
Формула (10) задает электрическое поле системы зарядов без учета граничных условий. Как правило, пучок движется в некотором ограниченном пространстве, что требует дополнить формулу (10) граничными условиями. Введение граничных условий приводит к дополнительному усложнению вычислительного процесса, так как в этом случае помимо вычислений по (10) приходится решать граничное интегральное уравнение (ГИУ) или краевую задачу для уравнения Лапласа.
Проиллюстрируем постановку такой краевой задачи. Пусть задана область Q, в которой ищется распределение поля E = —Vu, тогда можно записать систему
Дм (s) = —р (s), s £ Q, u\r = фо,
значение потенциала фо на границе Г области Q считается известным. Обычно это условие выполнено для широкого круга задач. Далее представим потенциал u в виде суперпозиции потенциалов up, которые находятся по формуле метода PP, и некоторого добавочного потенциала Su:
u = up + Su.
В этом случае справедливы следующие постановки задач:
Дup = р, up\r = фp;
ДSм (s) = 0, s £ Q,
Su\r = фо — фр.
Первую задачу решать не требуется, так как ее решение будет найдено из метода РР. Единственное, что потребуется, так это определить потенциал на границе Г. А вот вторую краевую задачу надо будет решать в области П разностными методами или воспользоваться ГИУ. Метод ГИУ основан на использовании формулы Грина и заключается в решении на границе Г области П интегрального уравнения вида
и(М) = £/и (р) - ± J ^^-и (р) dsp+± J
где тмр есть расстояние между точками M и p, а np - вектор внешней нормали к поверхности Г.
Функция потенциала на границе Г считается заданной, и требуется вычислить распределение поверхностной плотности заряда <т = ^. Зная плотность заряда на поверхности области, можно через интеграл рассчитать соответствующий ей потенциал.
Видно, что приведенный подход обладает достаточно большой вычислительной емкостью, чтобы подумать о существовании других методов, основанных на более «легковесных» алгоритмах. Тем более, что указанную задачу о нахождении собственного поля Es требуется решать на каждом шаге интегрирования по времени.
2г. Метод частиц-в-ячейках (Particle-In-Cell method, PIC). Этот метод обладает рядом преимуществ в сравнении с описанным выше методом PP. В данном подходе с целью понижения вычислительных затрат по определению поля Es используется решение краевой задачи для уравнений Максвелла. Идея состоит в том, что число узлов разностной! сетки, используемой при решении краевой задачи, существенно ниже, чем число O (N2). А численный метод решения краевой задачи с помощью быстрого преобразования Фурье (БПФ) [3] достаточно эффективен. Благодаря этим фактам удается получить выигрыш в скорости.
Опишем подробнее метод PIC. Пусть имеется некое распределение макрочастиц в пространстве V. Далее выделяется некая подобласть Q, содержащая рассматриваемый пучок и имеющая границу Г, на которой известны краевые условия типа Дирихле или Неймана для электростатического потенциала у. В этом случае из уравнений Максвелла (6) легко получить следующую краевую задачу:
А<р(р) = -£&, реп
дер дп
Фп, rD и rN = Г, (11)
где Е8 = -Ур.
Такая задача решается на каждом шаге интегрирования по времени уравнений движения частиц. Для решения задачи (11) используются разностные методы, т. е. такие величины как потенциал, плотность заряда, вектор электрического поля определены в узлах некой сетки, накинутой на область П. Эту сетку обычно называют эйлеровой. А положение макрочастиц, на которые действуют указанные поля, задается так называемыми лагранжевыми координатами.
В каждый момент времени есть некоторое распределение макрочастиц. Для него можно построить функцию плотности заряда, заданную в узлах, указанной разностной сетки. Здесь существует масса тонких моментов, касающихся способа «раздачи» плотности заряда в узлы сетки. Алгоритм представления облака плотности заряда макрочастицы зависит от параметров облака, и от его выбора будут меняться параметры функции плотности и соответственно самого электрического поля.
JV
Итак, для численного решения такой задачи на область Q накидывается разностная (эйлерова) сетка. Далее на полученной сетке находится функция плотности заряда. В результате искомое поле Es будет определено в узлах сетки.
В целом процесс математического моделирования поставленной задачи требует большого объема вычислительного времени и анализа получаемых при этом данных. Потому и встал вопрос об ускорении процесса моделирования с использованием массивно параллельных вычислений на графических процессорах GPU. Начиная с 2006 г. корпорацией NVIDIA развивается программная среда CUDA для массивно параллельных вычислений на GPU. Вычисления на GPU имеют отличительную особенность от распространенных вычислений на кластерах, состоящих из CPU. Графический процессор не создан для замены центрального процессора, он является, по сути, математическим сопроцессором. Особенностью этого GPU является то, что он содержит не 4 и не 8 ядер, а сотни, на уровне 240-512 ядер на одной видеокарте, которая может быть расположена в обычном персональном компьютере. Для примера, в один компьютер можно поставить до 4 таких видеокарт. Так уже существуют и серверные варианты таких GPU, например модель Tesla C1070 содержит 960 ядер. Из таких плат можно составить серверную стойку и, таким образом, получить огромный прирост в производительности по сравнению с используемыми CPU кластерными решениями. Вторым преимуществом рассматриваемого подхода является быстрый обмен данными между ядрами GPU процессора, так как все они локализованы на одном чипе, в отличие от кластерных CPU решений, на которых порой существенная часть времени теряется на обмене данными между процессорами.
Чтобы использовать ресурсы параллельных суперкомпьютеров, необходима разработка алгоритма распараллеливания физической задачи.
В качестве программы для расчетов была использована написанная ранее программа CBDA (Cyclotron Beam Dynamics Analysis) [4]. Она ориентирована на расчет циклотронов. Однако, с точки зрения математического моделирования, линия инжекции в циклотроне схожа с линией инжекции для синхротрона. В работе [5] было произведено переложение основного вычислительного алгоритма на GPU. И как следствие, скорость вычислений на Tesla C1060 выросла в 60-90 раз в зависимости от параметров модели. Такой коэффициент ускорения открывает принципиальные возможности для детальной оптимизации и моделирования синхротрона, особенно задач с большим пространственным зарядом.
3. Результаты математического моделирования поставленной задачи. Данная задача является частью проекта медицинского синхротрона ОИЯИ [6]. Компьютерное моделирование синхротрона осуществлено с помощью программ CELLY [7] и BeamLiner [8]. В нашей задаче требовалось произвести моделирование линии транспортировки пучка от источника ионов углерода 12С6+ с интенсивностью от 25 до 100 мА до RFQ структуры. Согласование параметров входа-выхода линии транспортировки без пространственного заряда проведено с помощью программы [9], математические алгоритмы для которой изложены в работе [10]. Параметры предложенной линии инжекции приведены на рис. 1.
На рис. 2 показаны распределения магнитного поля вдоль осевой линии канала транспортировки для каждого из обоих фокусирующих соленоидов. Суммарное магнитное поле в канале определялось как суперпозиция указанных полей.
2, мм
Соленоид 1
2сеп,ег = 395 ММ
Я тт = 60 мм
*тах = ЮОММ ^се„,ег = 2.3кГс
Рис. 1. Общий вид линии инжекции
В, кГс 2.5 |-
-200 0 200 400 600 800 1000 1200 1400 1600 1800
Z, мм
Рис. 2. Распределение магнитного поля вдоль осевой линии канала транспортировки (Тт = 73.2 кА)
Проекции начального эмиттанса пучка иллюстрирует рис. 3. Потенциал инжекции составляет 80 кВ. Необходимо было произвести расчеты для токов 0, 25, 50 и 100 мА. В силу большой величины тока из источника следует производить учет эффекта пространственного заряда, что, в свою очередь, требует моделирования большого числа частиц.
Результаты моделирования представлены на рис. 4. Здесь показано движение пучка от источника ^ = 0 мм) до входа в RFQ ^ = 1545 мм) для различных токов из источника: 0, 25, 50 и 100 мА.
Л, мрад
2оу = 28.972 мм >ру = 22.448 мрад = -2.12
Ру = 3.026
мм мрад
= 277.424л • мм • мрад
0 200 -40
я град. -20 0 20 40
X, мм
Рис. 3. Начальный эмиттанс из источника (см. рис. 1)
Рис. 4- Прохождение пучка через линию инжекции Ток пучка: а - 0 мА, б -25 мА, в - 50 мА, г - 100 мА.
1600
2, мм
Приведем получившиеся эмиттансы на входе в RFQ в зависимости от величины тока источника:
4. Заключение. Данная работа выполнена в рамках проекта ЛФВЭ ОИЯИ (г. Дубна) «Разработка и создание узлов ускорителя для адронной терапии». Для линии ин-жекции приведена математическая постановка задачи с учетом эффекта пространственного заряда и описаны алгоритмы ее решения.
Проведено компьютерное моделирование линии транспортировки пучка от источника ионов углерода с интенсивностью от 25 до 100 мА до RFQ структуры. Моделирование производилось программой CBDA с поддержкой вычислений на GPU. На основании моделирования определена пропускная способность линии инжекции.
Литература
1. Хокни Р., Иствуд Дж. Численное моделирование методом частиц / пер. с англ. А. С. Алпатова, А. Н. Полюдова; под ред. Р. З. Сагдеева, В. И. Шевченко. М.: Мир, 1987. 638 с. (Hockney Roger W., Eastwood Jeames W. Computer simulation using particles.)
2. CUDA Technology // URL: http://www.nvidia.com, http://www.nvidia.ru.
3. Рошаль А. С. Быстрое преобразование Фурье в вычислительной физике // Изв. вузов. Радиофизика. 1976. Т. XIX, № 10. С. 1425-1454.
4. Perepelkin E., Vorozhtsov S. CBDA — Cyclotron Beam Dynamics Analysis code // Proc. the XXI Russian Accelerator Conference (RuPAC2008). Zvenigorod, Russia, September 28 — October 3. Zvenigorod, Russia, 2008. P. 40—42.
5. Перепелкин Е., Смирнов В., Ворожцов С. Использование технологии NVIDIA CUDA при моделировании динамики пучка в ускорителях заряженных частиц // Вестн. Рос. ун-та дружбы народов. Сер. Математика. Информатика. Физика. 2010. № 1. С. 76—82.
6. Юдин И. П., Панасик В. А. Управление дополнительным каналом транспортировки пучка углерода при выводе из ускорителя Нуклотрон для медицинских целей // Устойчивость и процессы управления. Всерос. конференция, посвященная 80-летию со дня рождения В. И. Зубова. Санкт-Петербург, 1—2 июля 2010 г. СПб.: ВВМ, 2010. С. 134.
7. Юдин И. П. CELLY — программа расчета синхротронов // Программы расчета и моделирования для ускорительной техники. М.: ИТЭФ, 1992. С. 17—18.
8. Yudin I. P., Trofimov A. V. BEAMLINER — an object oriented beam line modeling C++ Code // Proc. PAC'99 (Particle Accelerator Conference). March 27—31. 1999. Brookhaven, USA, 1999. P. 1121.
9. Yudin I. P. Program for calculations of aberrations in solenoidal lens // Intern. Conference on Ion Sources (ICIS'03): abstracts. Dubna, Sept. 8—13, 2003. Dubna: JINR, 2003. P. 112.
10. Юдин И. П. Получение решений нелинейных уравнений математической физики с помощью метода функции влияния // Изв. Саратовск. гос. ун-та. Сер. Физика. 2010. № 1. С. 25—30.
Статья рекомендована к печати проф. Д. А. Овсянниковым. Статья принята к печати 26 апреля 2012 г.
Величина тока, мA
Радиус, мм ............
ar ..........................
вг, мм/мрад ..........
£r, ж- мм • мрад .....
0 25 50 100
5 20 38 50
0 1.5 1.6 1.3
0.13 1.4 3.2 5
276 330 410 580
