Для решения ряда задач физики плазмы (например, моделирование осесимметричных магнитных ловушек) целесообразно использовать цилиндрическую систему координат. Однако, при построении алгоритмов на основе метода частиц в ячейках, в цилиндрической геометрии возникает ряд особенностей, существенно снижающих эффективность параллельных вычислений. В данной работе на примере задачи моделирования динамики плазмы в осесимметричной ловушке-мишени рассматривается несколько алгоритмов, позволяющих устранить особенности, возникающие при использовании цилиндрической системы координат. С помощью предложенного подхода удалось достичь масштабируемости задачи до нескольких тысяч вычислительных ядер и рассчитать параметры мишенной плазмы во всей области ловушки. Ключевые слова: физика плазмы, метод частиц в ячейках, высокопроизводительные вычисления.
Особенности использования цилиндрической геометрии при решении задач физики плазмы методом частиц в ячейках* 2. Математическая модель 2.1. Плазменная ловушка-мишень
Рассмотрим построение математической модели для описания динамики плазмы в ловушке-мишени. Мишенная плазма представляет собой цилиндрический столб плазмы, удерживаемый в магнитной ловушке магнитными полями достаточной длины и апертуры для прохождения нейтрализуемого пучка отрицательных ионов. Удержание плазмы магнитным полем представляет собой непростую задачу. Требуется создать магнитное поле, способное удерживать мишенную плазму, и геометрия которого способна подавлять неустойчивости плазмы. Особо необходимо хорошее удержание плазмы от истечения в торцевые отверстия-пробки, иначе вытекающая через них плазма вынесет большое количество газа, и потребуется чрезвычайно производительная система откачки. В осесимметричной магнитной ловушке в качестве торцевых пробок используются инверсные пробки, в которых создается продольное магнитное поле с противоположным направлением относительно поля во всей ловушке между пробками. С помощью мультипольных магнитных стенок, сформированных последовательностью кольцевых магнитов с переменной намагниченностью, формируется поле, препятствующее радиальным потерям плазмы.
Модель ловушки представляет собой цилиндрическую систему длиной 1355 мм с вакуумной камерой внутренним диаметром 199 мм. Ловушка состоит из центрального катодного блока и симметричных относительно него цилиндрических секций и торцевых крышек. Схема ловушки приведена на рисунке 1.
r rotE = − ,
r c ∂t r divE = 4πρ = 4π ∑qα ∫ fα dv ,
v α divB = 0 .
Рис. 1. Схема плазменной ловушки-мишени. (1) – широкоапертурные отверстия для прохождения нейтрализуемого пучка, (2) – катодный блок, (3) – один из шести катодов, (4) – постоянные магниты, (5) – магнитный экран.
В ловушке с магнитным полем происходит запирание электронов, поэтому важно в расчетах учитывать фазовое распределение, область потерь и область удержания. Математическое моделирование необходимо для оценки и минимизации потерь плазмы в проходные отверстия в торцах (через инверсные пробки), а также через мультипольные магнитные стенки ловушки на ее вакуумную камеру. При этом наиболее полно динамика плазменных электронов может быть описана уравнением Больцмана для функций распределения ∂fα + vr ∂frα + qα (Er + 1 [vr, Br]) ∂frα = St{ fα } , (1) ∂t ∂r c ∂p и системой уравнениrй Максвелла с самосогласовrанными электромагнитными полями [2]: r 4π r 1 ∂E 4π r r 1 ∂E rotB = j + = ∑ qα ∫ fα vd v + , (2) c r c ∂t с α c ∂t 1 ∂B (5) r r Здесь индексом α обозначается сорт частиц (ионы и электроны), p = mγ v - релятивистский импульс частицы, St{ fα } - функция, описывающая следующие физические процессы: • • • ионизация атомов водорода катодными электронами H + e- → H+ + 2e-, ионизация молекул водорода катодными электронами H2 + e- → H2+ + 2e-, ионизация ионами H+ атомов водорода H2 + H+ → H3+ + H 2.2. Решение основных уравнений
Для решения системы (1-5) используется комбинация метода частиц в ячейках [3] и метода Монте-Карло. Чтобы описать физические эффекты, происходящие во всей ловушке с высокой точностью, необходимо использовать не менее 10^9 модельных частиц и не менее 10^6 узлов сетки, что требует применения высокопроизводительной техники. Такие требования к модели вызваны сложной конфигурацией магнитного поля, многокомпонентностью плазмы, а также необходимостью учитывать ионизацию газа и истечение плазмы из ловушки.
Решение уравнения Больцмана методом расщепления сводится к решению уравнения Власова и корректировке траекторий частиц с учетом физических процессов с использованием метода Монте-Карло [4].
Движение частиц происходит согласно характеристикам уравнения Власова по leap-frog схеме,rобеспечивающей второrй порядок точности: dp r 1 r r dr r
= qα (E + [v, B]) , = v . (6) dt c dt prm+1/ 2 − prim−1/ 2 i 1 vrim+1/ 2 + vrim−1/ 2 , Brim , c = qα Erim + rim+1 − rim = vim+1/ 2 . (7) τ 2 τ Здесь τ - шаг по времени; верхний индекс указывает на момент времени, в который вычисляется искомая функция, нижний индекс i указывает на номер частицы, для которой производятся вычисления.
Для того чтобы избежать особенностей у оси, используется локальное преобразования из декартовых координат в цилиндрические – схема Бориса [5].
Для моделирования процесса ионизации используется вероятность столкновения заряженной частицы с нейтральными атомами. Если столкновение произошло, в расчетную область добавляются частицы с необходимыми параметрами.
Вероятность столкновения для j-ой частицы, имеющей скорость vj, за время ∆t может быть рассчитана по формуле [6] v
Pj = 1 − exp(−∆s jσ (ε j )n j (rj )), (
Решение уравнений Максвелла (2-3) осуществляется по схеме Лангдона-Лазински [7]: r r Bm+1/ 2 − Bm−1/ 2
r = −c rothEm , r τ r Em+1 − E m τ r divh Evm = 4πρ m , divh Bm+1/ 2 = 0 .
r r
= −4π j m+1/ 2 + c rothBm+1/ 2 ,
(9)
(
j r r Здесь q j - заряд частицы с номером j; функция R(r , rj (t)) (функция ядра) характеризует форму, размер частицы и распределение в ней заряда.
Как показано в [8], при таком способе вычисления плотности тока и плотности заряда
уравнения (
+ divh j m+1/ 2 = 0.
τ
Поэтому в общем случае необходимо скорректировать значение плотности тока для
выполнения уравнения (15). Этот метод описан в [5].
r
На первом этапе определим электрическое поле E по схеме (
r ∆hδ Φ = divhE* − ρ . r После этого можно получить окончательное значение E и приступить к нахождению магнитного поля по схеме (9).
Такой метод является более ресурсоемким, поскольку требует решение уравнения Пуассона, кроме того, существенно усложняет построение параллельного алгоритма на основе декомпозиции области.
Для того чтобы избежать решения уравнения Пуассона и корректно находить сеточные значения напряженности электрических и магнитных полей, необходимо изначально потребовать выполнения уравненrия неразрывности (15), тогда уравнение (11) выполнено в любой момент времени, если divhE0 = ρ 0 . В работах [9-11] для некоторых функций R вычисления плотности тока и плотности заряда были согласованы в соответствии с (15). Однако, основным условием такого согласования является симметрия ядра
R((x, y, z), (x j (t), y j (t), z j (t))) = S(x, x j (t)) * S ( y, y j (t)) * S (z, z j (t)) .
В цилиндрической системе координат такой симметрии нет и функция R в узле с номером (i, l, k) имеет вид: (13) (14) (15) (16) 1 * ri2+1 − rj2 * hϕ − |ϕ l −ϕ j | * hz − | zk − z j | , Vi ri2+1 − ri2 hϕ hz ri < rj < ri+1,|ϕ l −ϕ j |< hϕ ,| zk − z j |< hz .
R((ri ,ϕ l , zk ), (rj,ϕ j, z j)) = 1 * rrii22−−11 −− rrij22 * hϕ − |ϕ l −ϕ j | * hz − | zk − z j | , Vi hϕ hz ri−1 < rj < ri ,|ϕ l −ϕ j |< hϕ ,| zk − z j |< hz . 0, иначе где Vi = 2πrihrhϕ hz - объем ячейки. Мы не рассматриваем здесь случай r=0, поскольку он не влияет на дальнейшее построение параллельного алгоритма и может быть рассмотрен независимо. 1 tm+1
∫ ( − vrm+1/ 2 d (SR(ri , rj (t)))* SP(ϕ l ,ϕ j (t)) * SZ (zk , z j (t)))dt + τVi tm dr 1 tm+1
∫ ( − vϕm+1/ 2 d (SP(ϕ l ,ϕ j (t)))* SR(ri , rj (t)) * SZ (zk , z j (t)))dt + τVi tm dϕ 1 tm+1
∫ ( − vzm+1/ 2 d (SZ (zk , z j (t)))* SR(ri , rj (t)) * SP(ϕ l ,ϕ j (t)) )dt τVi tm dz
2 ∆f = f m+1 − f m , δf = Пусть частица переместилась из точки (r m ,ϕ m , zm ) в точку (r m+1,ϕ m+1, z m+1 ). Введем следующие операторы: f m+1 + f m SR(ri , rj (t)) * SP(ϕ l ,ϕ j (t)) * SZ (zk , z j (t))dt = (17) тая,
В данной работе мы рассмотрим вычисление функции плотности тока через прямое решение уравнение неразрывности. Поскольку в исходной постановке задачи необходимо учитывать все три компоненты скорости, то и уравнение неразрывности необходимо рассматривать в трехмерном пространстве, с последующим устранением зависимости отϕ .
Рассмотрим функцию ядра частицы (16) как комбинацию трех одномерных ядер 1 R((r,ϕ , z), (rj (t),ϕ j (t), z j (t))) = SR(r, rj (t)) * SP(ϕ ,ϕ j (t)) * SZ (z, z j (t))
V И запишем уравнение неразрывности в узле (i, l, k) для промежутка времени [t m ,t m+1], счичто частица при перемещении не покидает окрестность и движется с постоянной скоростью и выполним переход от [t m ,t m+1] к − 1 , 1 .
2 2 Траектория частицы тогда описывается следующим образом: r(t) = δr + ∆rt ϕ (t) = δϕ + ∆ϕ t z(t) = δz + ∆zt Подставив эти выражения в уравнение (17), получим 1 1∫/2 vrm+1/2 2r(t) SP(ϕl ,ϕ j (t))SZ(zk , zj (t))dt + Vi −1/2 ri2+1 − ri2 + =
V1i −11∫//22vϕm+1/2 h1ϕ SR(ri ,rj (t))SZ(zk , zj (t))dt + V1i −11∫//22vzm+1/2 h1z SP(ri ,rj (t))SR(zk , zj (t))dt = 2*vrm+1/2
*δr *(ϕl+1 −δϕ)(zk+1 −δz) + Vi *(ri2+1 − ri2)hϕ hz ∆r∆ϕ(zk+1 −δz) + ∆z∆ϕδr + ∆r∆z(ϕl+1 −δϕ)
+ 12
(18) + + JZim,k++11//22 = Vi *v(rzmi2++11/−2 ri2 ) * (ri2+1 −δr *δr) + ∆1r2∆r где Vi = 2πrihrhz - объем ячейки.
В случае, если частица перемещается между ячейками, то можно использовать алгоритмы распределения плотности тока между ячейками, описанные в [8,12]
Заметим, что уравнение (17) справедливо для любого ядра конечного размера, поэтому все выкладки, полученные выше, можно обобщить на все типы конечного размера ядер и любую систему координат.
Таким образом, решение уравнений (1-5) разбивается на два этапа – на первом,
лагранжевом этапе происходит движение частиц по схеме (7), с учетом столкновений (
Прямое вычисление плотности тока по формуле (20) требует решения уравнения Пуассона только в начальный момент времени и в дальнейшем все сеточные величины зависят только от значений в соседних узлах сетки на предыдущем временном шаге. Это дает возможность выполнить декомпозицию области, сводя все коммуникации только между границами подобластей, и достичь высокой масштабируемости параллельных алгоритмов.
Таким образом, можно построить следующий алгоритм параллельных вычислений, основанный на эйлерово-лагранжевой декомпозиции. Область решения разбивается вдоль координаты Z на несколько подобластей. Каждой подобласти выделяется группа процессоров и частицы подобласти равномерно распределяются между процессорами этой группы независимо от координаты. На каждом временном шаге внутри подобласти независимо вычисляются траектории частиц, определяется плотность тока на сетке. Внутри группы все сеточные значения суммируются. Граничные элементы сетки, а также частицы, покинувшие подобласть, пересылаются между группами. На основе этого подхода, авторами статьи разработан параллельный алгоритм с учетом статической балансировки вычислительной нагрузки, позволяющий добиться высокой масштабируемости до нескольких тысяч вычислительных ядер [13]
Как показано в работе [13], при расчете динамики 5*109 частиц на сетке 512х1024 ускорение при использовании 128 ядер по сравнению с 1024 составляет 7,71. 128 256 512 1024 На данной задаче также была проверена масштабируемость полученного алгоритма относительно 128 процессорных ядер. Результаты вычислительных экспериментов приведены в таблице 1.
Параметры расчетов: Число частиц 2,5*109, размер сетки 512х1024.
Расчеты проводились с использованием суперкомпьютеров – «Ломоносов» (Московский государственный университет), NKS-30T (Сибирский суперкомпьютерный центр СО РАН), кластера Новосибирского государственного университета.
Таблица 1. Масштабируемость алгоритма на различных архитектурах Количество использу- Время расчета одного шага, с. / ускорение относительно 128 ядер емых процессорных ядер Ломоносов МГУ, Intel Xeon X5570 2932 Mhz 5,60 / 1,00 2,96 / 1,89 1,51 / 3,70 0,76 / 7,37
Кластер НГУ, Intel Xeon X5670 2932 Mhz 9,24 / 1,00 4,92 / 1,88 2,46 / 3,76
НКС-30Т, Intel Xeon E5540 2530
Mhz 12,80 / 1,00 6,70 / 1,91 3,39 / 3,78 Как видно из таблицы, несмотря на то, что масштабируемость алгоритма на суперкомпьютере Ломоносов чуть хуже, время расчета одного шага существенно меньше из-за более быстрой работы с памятью. В целом, полученные результаты позволяют говорить о достижении высокой масштабируемости параллельного алгоритма на различных суперкомпьютерах в пределах одной тысячи процессорных ядер. В дальнейшем планируется увеличить число используемых процессоров до десятка тысяч по аналогии с [13] и использовать ускорители Intel Xeon Phi. 3. Результаты моделирования
Для адекватного воспроизведения магнитного поля ловушки требуется шаг пространственной сетки h не более 0,1 мм. Таким образом, при размере области 1355 мм х 99.6 мм, выбрана сетка 13550х996 узлов. В расчетах использовалось до 5*109 модельных частиц. Была рассчитана плотность ионов H+ во всей ловушке (Рис. 2)
Рис. 2. Плотность ионов H+ в ловушке, n_i = 5,72*1012 см-3. Момент времени 10-5 с. Из рисунка 2 видно, что ионы распространяются до самой оси ловушки, обеспечивая высокую плотностью плазмы. При этом значение плотности ионов по оси может отличаться в несколько раз, но интеграл плотности по осевой координате должен соответствовать требуемому значению. Также важно, чтобы в области пролета пучка (до 5 см от оси ловушки) радиальная плотность частиц была достаточно высокой. Ионы плазмы должны быть однородно распределены по радиусу для максимально эффективного использования пучка. Для проверки этого условия была рассчитана плотность ионов H+ вдоль всей длины ловушки на разном удалении от центра ловушки (Рис. 3.) Рис. 3. Плотность ионов H+ в области пролета пучка для различных значений радиуса, n_i = 5,72*1012 см-3. Момент времени 10-5 с.
Из рисунка 3 видно, что интегральные значения плотности по оси при удалении от оси существенно не изменяются. Таким образом, эффективно нейтрализуется весь пучок отрицательных ионов, проходящий через ловушку. Кроме того, важным моментом является резкое падение плотности частиц у торцов ловушки - это соответствует высокой степени удержания частиц.
Также для оценки эффективности удержания плазмы было рассчитано распределение частиц в области магнитных пробок ловушки.
Радиальное распределение плотности ионов в области магнитных пробок (Z=1296 мм) приведено на рис 4.
Рис. 4. Радиальная плотность ионов H+, H2+, H3+, n_i = 5,72*1012 см-3.
Как видно из рисунка, плотность ионов быстро падает по радиусу, при этом ионы H2+ и H3+ более размазаны по радиусу. Тем не менее, при радиусе 25 мм концентрация ионов различных типов выравнивается и стремительно падает. Это означает, что наружу проникают только частицы с повышенной поперечной энергией. Аналогичный результат был получен экспериментально в работе [14].
Результаты расчетов плотности плазмы на оси в области магнитной пробки представлены на рис. 5. Наблюдается ступенчатое падение плотности всех ионных компонент плазмы вдоль оси в области инверсной магнитной пробки. Этот результат также полностью соответствует наблюдениям работы [14]. Рис. 5. Плотность ионов H+, H2+, H3+ на оси, n_i = 5,72*1012 см-3и распределение Z-компоненты магнитного поля.
Достигнутая в расчетах средняя плотность плазмы составила до 2*1013 см-3. Именно такой плотности предполагается достичь на создаваемой ловушке. 4. Заключение
Разработана математическая модель ловушки-мишени для высокоэффективной нейтрализации мощных пучков отрицательных ионов. Модель построена на основе комбинации метода частиц в ячейках и метода Монте-Карло и описывает удержание столба плазмы в вакуумной камере магнитным полем сложной геометрии. Для расчета плотности тока в цилиндрической системе координат на основе уравнения неразрывности были получены прямые формулы, не требующие корректировки электрического поля с помощью уравнения Пуассона. Это позволяет избавиться от глобальной зависимости по данным и построить масштабируемый до нескольких тысяч вычислительных ядер параллельный алгоритм.
Использование современных суперкомпьютеров позволило рассчитать динамику плазмы, определить потоки плазмы в торцевые отверстия ловушки. Результаты расчетов показали, что образующаяся плазма имеет достаточную степень однородности и заполняет практически всю длину ловушки.
С помощью численных экспериментов показано, что магнитная система со слабым продольным полем и инверсными пробками в торцевых отверстиях в магнитном поле позволяет добиться достаточно малого потока плазмы из ловушки.
Также показано формирование на оси ловушки достаточно широкого (около 50 мм в диаметре) столба плазмы, имеющего плотность порядка 1013 см-3
Таким образом, на основе построенной математической модели и эффективных вычислительных алгоритмов еще до проведения лабораторных экспериментов удалось сделать важные оценки характеристик плазменной ловушки. В дальнейшем планируется взаимное уточнение результатов моделирования и эксперимента, оптимизация магнитной системы ловушки для достижения ее высоких параметров. Литература 7. Langdon A.B, Lasinski B.F. Electromagnetic and relativistic plasma simulation models // Meth.
Comput. Phys. 1976, Vol. 16. P. 327-366. 11. Regine Barthelme. Conservation de la charge dans les codes PIC. // Comptes rendus – Mathematique. 2005, Vol. 341, Issue 11. P 689-694 12. T. Umeda, Y. Omura, T. Tominaga, H. Matsumoto, A new charge conservation method in electromagnetic particle-in-cell simulations, Computer Physics Communication 2003, 156, Issue 1, P. 73-85. 14. Dimov, G. I., Emelev, I. S. Experiments to study the confinement of a target plasma in a magnetic trap with inverse plugs and circular multipole walls. // Technical Physics, 2014, Vol. 59, Issue 2, P. 181-189. The features of using the cylindrical geometry for solution of the plasma physics problems with the particle-in-cell method
E.A. Berendeev, M.A. Boronina, V.A. Vshivkov, A.A. Efimova Institute of Computational Mathematics and Mathematical Geophysics SB RAS A big number of problems in plasma physics (for example, for axially symmetric magnet traps) reasonably requires the usage of the cylindrical coordinate system. However, in the algorithms based on the particle-in-cell method, there are many peculiarities of the cylindric geometry, which significantly reduce the effectiveness of the parallel computations. In the present work, few algorithms are considered with the aim to reduce the problems of the usage of the cylindrical coordinate system. As an example, the modeling of the plasma dynamics in the axially symmetric trap is considered. The new approach allowed to reach the scalability for few thousands of processor cores. The plasma parameters in the whole domain were found. 3. C. K. Birdsall, A. B. Langdon. Plasma Physics via Computer Simulation. McGraw-Hill, 1985. 4. 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. P. 65-83. 5. Boris J.P. Relativistic plasma simulation - optimization of a hybrid code // Fourth Conference on nu-merical Simulation of Plasmas. Washington, 1970. P. 3-67. 6. V. Vahedi, M. Surendra. A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges // Comput. Phys. Comm., 87 (1995), P. 179-198. 7. Langdon A.B, Lasinski B.F. Electromagnetic and relativistic plasma simulation models // Meth.
Comput. Phys. 1976, Vol. 16. P. 327-366. 9. T. Zh. Esirkepov, Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor, Computer Physics Communications 2001. Vol. 135. P. 144-153. 11. Regine Barthelme. Conservation de la charge dans les codes PIC. // Comptes rendus – Mathematique. 2005, Vol. 341, Issue 11. P 689-694 12. T. Umeda, Y. Omura, T. Tominaga, H. Matsumoto, A new charge conservation method in electromagnetic particle-in-cell simulations, Computer Physics Communication 2003, 156, Issue 1, P. 73-85. 13. Berendeev E.A., Boronina M.A., Korneev V.D. Parallel'nyy algoritm resheniya zadach di-namiki zaryazhennykh chastits s uchyetom balansirovki vychislitel'noy nagruzki. (Parallel algorithm for charge particles dynamics with load-balance calculations)// Vestnik YuUrGU. Seriya “Vychislit