CC BY
74
Вычислительные технологии Том 6, № 2, 2001
АЛГОРИТМЫ РЕШЕНИЯ ЗАДАЧИ ВЗАИМОДЕЙСТВИЯ ЛАЗЕРНОГО ИМПУЛЬСА С ПЛАЗМОЙ *
В. А. Вшивков
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected] К. В. Вшивков Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия
Г. И. ДудниковА
Институт вычислительных технологий СО РАН, Новосибирск, Россия
Two-dimensional and three-dimensional algorithms for the solution of the problem of laser pulse interaction with plasma are considered. The full kinetic model for plasma description is used. The problem is solved by the particle-in-cell (PIC) method. Several approaches to parallel realization of the algorithm with various degrees of processor units loading are investigated. The results of the modelling are presented.
Введение
Последние достижения в развитии лазерной техники позволили получить сверхкороткие, длительностью в несколько фемтосекунд (1 ф. с = 10-15 с), импульсы электромагнитного излучения с интенсивностью 1021 Вт/см2 [1]. В связи с этим встала проблема изучения взаимодействия лазерных импульсов с веществом, в том числе с плазмой. Одним из методов исследования этой проблемы является численное моделирование [2-5]. Вычислительная сложность данной задачи связана с наличием нескольких существенно разных временных и пространственных масштабов. В пространстве — это длина волны лазерного импульса, которая должна быть описана достаточно хорошо (на длину волны надо иметь до 10 ячеек), и пространственный размер задачи, который достигает нескольких десятков длин волн лазерного импульса. Во времени — это существенно разные скорости движения ионов и электронов. Кроме того, многие процессы требуют трехмерного моделирования. Только в последнее время с появлением суперкомпьютеров появилась возможность решать задачи, которые характеризуются большим объемом данных и вычислений. В настоящей статье представлены несколько подходов к параллельной реализации алгоритма решения задачи о взаимодействии короткого лазерного импульса с плазмой. Рассматриваемые подходы
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, гранты №99-01-00512, №99-07-90418 и Федеральной целевой программы "Интеграция", проект №274.
( В. А. Вшивков, К. В. Вшивков, Г. И. Дудникова, 2001.
характеризуются различной степенью загрузки процессорных элементов (ПЭ). Разработанная программа настраивается на доступное число ПЭ.
При взаимодействии лазерного импульса с тонким слоем плотной плазмы формируются импульсы с резким передним и/или задним фронтом [2]. Такие импульсы создают так называемую "кильватерную волну", в которой заряженные частицы ускоряются до релятивистских скоростей. Эти высоко энергетичные частицы могут в дальнейшем использоваться для решения проблемы управляемого термоядерного синтеза.
1. Постановка задачи
Рассмотрена задача в следующей постановке. В области, имеющей форму параллелепипеда, находится фольга, моделируемая как тонкий слой плазмы, состоящей из электронов и ионов одного типа. Через поверхность параллелепипеда входит лазерный импульс в виде пакета электромагнитных волн различной амплитуды и поляризации. Импульс взаимодействует с фольгой, частично отражаясь и проходя через нее. Режимы проникновения и отражения лазерного импульса зависят от его характеристик и плотности плазмы.
Область характеризуется линейными размерами, лазерный импульс — амплитудой волн, линейными размерами и положением в области. Плазма представлена набором достаточно большого числа модельных частиц с заданной плотностью (рис. 1).
Рассматриваемый физический процесс описывается уравнениями Власова для ионов и электронов, которые имеют вид
дДе , дДе
-- + + ^ье--
т
д г
д р
0, ¥г,е = %е[ Е + 1 ["V, В]
(1)
Уравнения (1) дополняются уравнениями Максвелла, связывающими между собой плотности заряда р и тока ^ напряженности электрического и магнитного полей:
„ 4ж. 1 дЕ В = — j + ,
с с дЬ
го1 Е = -
1 д В
С ~дЬ:
Е = 4пр,
в = о.
(2)
(3)
(4)
(5)
Входящие в уравнения Максвелла плотности тока и заряда вычисляются из функций распределения ионов и электронов по формулам
j = ^ Ф,е / /г,еУ (6)
Р = ^ ?г,е / Де
Здесь индексы % и е обозначают тип частиц — ионы и электроны соответственно; Де — функция распределения частиц; ^ — заряд; Е — напряженность электрического поля; В — напряженность магнитного поля; р = — релятивистский импульс частиц, 7 = (1 — V2/с2) 1/2, т - масса покоя иона или электрона.
Перейдем во всех уравнениях к безразмерным переменным, введя в качестве масштабов величины, указанные ниже.
Длина Длина волны лазерного импульса Л
Время Период волны лазерного импульса, Ь0 = 2п/^0,
где — частота электромагнитного импульса
Скорость Скорость света с
Напряженность электрического поля Е0 = тес^0/2пе
Плотность плазмы п0 = те^0/16п3е2
Уравнение Власова (1) в безразмерной форме примет вид
д/1,е д/1,е д/1,е п
дЬ от др (7)
Ее = — (Е + [V, В]), Ег = к (Е + [V, В]),
где к = те/т^ — отношение масс электронов и ионов; р = 7V, 7-1 = д/1 — V2. Уравнения Максвелла (2) - (5) преобразуются к виду
дЕ „ . дВ
ж = го1 В—^ ж = —10 Е,
Е = р, div В = 0.
Плотность заряда р и плотность тока j удовлетворяют уравнению неразрывности
жj = 0В начальный момент времени частицы, которые моделируют плазму, располагаются в области х0 < х < х0 + 8. Нами рассмотрены три основных варианта распределения частиц:
— частицы образуют тонкую фольгу (8 ~ Л);
— частицы представляют собой две последовательно расположенные тонкие фольги;
— частицы распределены равномерно во всей области моделирования (8 = хтах). Частицы имеют разные заряды (ионы — положительный, электроны — отрицательный); скорости частиц, компоненты напряженностей электрического и магнитного полей равны
нулю. Задана плотность плазмы. Лазерный импульс входит через левую границу области, известны его линейные размеры и амплитуда волн.
На левой границе оси x задаются условия входа импульса в область, а на правой — условия выхода. Условия на границах области в направлениях осей y и z берутся периодическими, т.е. f (t,x, 0, z) = f (t,x,ymax,z), f (t,x,y, 0) = f (t,x,y, zmax), где f — любая функция из уравнений; ymax, zmax — линейные размеры области по осям y, z.
2. Методы и алгоритмы решения
Физические задачи решаются с помощью сложной системы дифференциальных уравнений для функций распределения ионов и электронов по скоростям, которые являются функциями семи аргументов (времени, трех пространственных переменных и трех компонент скорости), и уравнений Максвелла для функций четырех аргументов (времени и трех пространственных переменных).
Решение уравнения Власова. С математической точки зрения уравнение Власова является уравнением гиперболического типа, которое может быть записано в эйлеровых или лагранжевых координатах. Использование эйлеровых координат приводит к конечно-разностным схемам, в которых значения сеточных функций в соседних узлах введенной расчетной сетки связаны друг с другом аппроксимацией производных.
В лагранжевых координатах в шестимерном пространстве также можно ввести расчетную сетку, в узлах которой заданы значения функций. В переменных Лагранжа координаты узлов связаны с движением плазмы, что обычно является большим недостатком при решении уравнений конечно-разностными методами, потому что в ходе вычислений ла-гранжева сетка деформируется. Это делает невозможной аппроксимацию пространственных производных.
Уравнение Власова в лагранжевых координатах имеет вид df/dt = 0, т.е. в нем отсутствуют производные по пространству и скоростям. Решение его (f = const) говорит о том, что значение функции распределения не меняется при движении вдоль характеристики уравнения и не зависит от значения функции распределения в других точках пространства.
На основе этого представления уравнения Власова в лагранжевых переменных разработан метод частиц в ячейках — наиболее универсальный и широко применяемый в настоящее время метод для решения задач физики бесстолкновительной плазмы [6-8]. В этом методе вся плазма моделируется набором отдельных частиц, каждая из которых характеризует движение многих физических частиц. Характеристики уравнения Власова (7) описывают траектории этих частиц.
Уравнения этих характеристик имеют вид:
для электронов для ионов
I = -(E + vB])
dp
dt
k(E + [v, B]),
dr
dt
В настоящей работе для уравнений (8) используется схема
m+1/2 m— 1 /2
Р i - P i
T
(9)
Количество выбранных частиц определяется необходимой точностью решения и имеющимися ресурсами ЭВМ (памятью и быстродействием). Таким образом, метод частиц сводит задачу решения уравнения Власова в шестимерном пространстве к задаче интегрирования уравнений движения для отдельных частиц. Другими словами, временная сложность решения задачи уменьшается от O(n6) при решении конечно-разностными методами, до O(n3) в методе частиц.
Моделирование методом частиц носит очевидный имитационный характер. Модельные частицы двигаются по тем же траекториям и с теми же скоростями, что и реальные ионы и электроны. Количество модельных частиц в малом объеме, как и в реальной плазме, характеризует плотность этих частиц и т. д. Поэтому свойства реальной и модельной плазм подобны. Однако некоторые характеристики плазмы существенно отличаются от характеристик в реальной плазме. Например, в модельной плазме длина свободного пробега много меньше, а частота столкновений много больше, чем в реальной [6]. В модельной плазме также значительно выше уровень тепловых шумов. Для ликвидации этих недостатков в методе частиц используются частицы конечного размера.
Метод частиц имеет еще одно важное технологическое достоинство, заключающееся в том, что движение частицы в каждой ячейке зависит только от значений напряжен-ностей электрического и магнитного полей в ней и не зависит от координат и скоростей других частиц. Это дает возможность независимо рассматривать такие ячейки вместе с находящимися в них частицами, что позволяет легко адаптировать программу, реализующую метод, для решения различных задач. Изменению подвергаются только геометрия пространства моделирования, начальные и граничные условия, решение уравнений Максвелла. Движение же частиц внутри ячеек происходит по одним и тем же законам для всех задач. Независимость движения частиц друг от друга допускает распараллеливание алгоритма решения уравнений движения частиц и вычисления плотности плазмы и средних скоростей частиц. Распараллеливание алгоритма дает возможность увеличить количество модельных частиц, а следовательно, повысить точность расчетов, что позволяет исследовать более тонкие физические процессы.
Сравним объем оперативной памяти, необходимый для решения уравнения Власова при помощи конечно-разностных схем и PIC-метода. Пусть мы используем разбиение шестимерного пространства по каждому из пространственных или скоростных направлений на 25, т.е. 32 ячеек. Тогда длина массива для хранения значений функции распределения составит 230 ^ 109 слов, что неприемлемо по ресурсам. Для PIC-метода рассмотрим такую же пространственную сетку, что дает 215 ячеек. Пусть в каждой ячейке в среднем находятся 25 частиц, каждая из которых характеризуется семью числами (зарядом или массой, тремя пространственными координатами и тремя скоростями). В результате получаем 7 • 220 ^ 7 • 106 слов, что примерно в 150 раз меньше, чем в первом случае. На больших сетках эффект становится еще заметнее.
Решение уравнений Максвелла. Уравнения Максвелла решаются в эйлеровых переменных. Необходимые для их решения плотности заряда и тока определяются по ско-
ростям и координатам отдельных частиц:
рМ) = 5]^^(*)), -¡М) = 5] ^V(*№,г-(*))• ¿=1 ¿=1
Здесь — заряд частицы с номером ^; функция Л (г, г^- (¿)) характеризует форму, размер частицы и распределение в ней заряда [6].
Для нахождения электрических и магнитных полей обычно используется схема, предложенная Лэнгдоном и Лазинским в 1976 г. [9], в которой поля определяются из разностных аналогов законов Фарадея (2) и Ампера (3):
gm+1/2 _ gm-1/2
Т
—c rothEm
Em+1 — Em
Т
= -4njm+1/2 + c roth Bm+1/2.
(10) (11)
В этой схеме значения компонент напряженностей электрического и магнитного полей вычисляются в узлах сеток, смещенных относительно друг друга по времени и пространству (рис. 2). Разностные операторы roth и divh на указанной сетке имеют вид
rothB
Bz,i,fc,l-1/2 — Bz,i,k-1,1-1/2 By,i,fc-1/2,l — By,i,fc-1/2,l-1
. hy hz
Вж,г-1/2,k,l — Вж,г-1/2,k,l-1 Bz,i,fc,1-1/2 — Bz,i-1,fc,l- 1/2
hz hx
By,i,fc-1/2,l — By,i-1,fc-1/2,l -1,l 1
divhB
hB =
BXt+1/2,fc,l — Вж,г-1/2,k,l By,i,fc+1/2,l — By,i,fc-1/2,l Bz,i,fc,l+1/2 — Bz,i,fc,l-1/2
hx
+
hy
+
h2
Эта схема имеет второй порядок аппроксимации по пространству и времени, и, что особенно важно, в ней автоматически удовлетворяются разностные аналоги уравнений (4) и (5) при отсутствии зарядов:
а^вт+1/2 = 0,
divhEm = 0.
(12)
В случае, когда плотность заряда не равна нулю, все зависит от согласования вычислений плотностей заряда и тока. Если они вычисляются так, что удовлетворяется разностный аналог уравнения неразрывности
Р-^ + divh jm+1/2 = 0, (13)
т
то закон Гаусса также выполняется. В случае невыполнения этого уравнения электрическое поле подправляется так, чтобы его разностный аналог выполнялся. Такая поправка предложена Борисом в 1970 г. [10]. Она требует решения уравнения Пуассона на каждом шаге по времени, что отнимает много времени счета. В 1987 г. Мардер [11] предложил добавить в правую часть уравнения (11) дополнительный член, который ведет к уменьшению невязки в уравнении (12):
Em+1 _ pm
-= -4njm+1/2 + c rothBm+1/2 + d gradF,
т
где
F = divhEm - 4npm
Этот метод более экономичный, в нем не требуется решения уравнения Пуассона. Но его недостаток заключается в том, что на каждом шаге погрешность в уравнении (12) полностью не устраняется, что постепенно может искажать физические процессы.
В 1992 г. Вилласенор и Бунеман [12] предложили метод вычисления плотности тока при помощи точного учета потоков плазмы через границы ячеек. Это позволило автоматически удовлетворить разностному уравнению неразрывности (13) и, следовательно, точно выполнить разностный закон Гаусса. Аналогичный метод использован в наших программах.
Покажем возможность построения такого алгоритма на одномерном примере. Для уравнения (13) рассмотрим дивергентную разностную схему вида
m+1 m -m+1/2 -m+1/2
Pi—1/2 - Pi—1/2 + J ' - Ji-1 = 0 (14)
т h
Здесь сеточная функция плотности {pi—1/2} определена в узлах — центрах ячеек, а потоки — на границах ячеек.
Такой выбор дает второй порядок аппроксимации по h и обеспечивает выполнение сеточного аналога формулы Остроградского — Гаусса, т.е. консервативность. В схеме (14) также используются сдвинутые временные сетки, что позволяет получить второй порядок аппроксимации по т. Подобная схема распространена в задачах физики плазмы, когда для уравнения движения частиц применяется схема с перешагиванием (leapfrog scheme). В данном случае для определения координат частиц используется схема
Xm+1 _ xm
XP_= vm+1/2
т
Поскольку схема (14) линейна относительно входящих в нее сеточных функций, достаточно добиться выполнения разностного уравнения неразрывности для каждой отдельной частицы. Тогда его справедливость для полных плотностей и потоков, получаемых суммированием по всем частицам, будет гарантирована.
Плотности заряда и тока определим формулами
Р^ = 5] 3РЛ(Жг-1/2 - О, (15)
(р)
зт+1/2 = £ ^+1/2я(хг - х;+1/2). (16)
(р)
Для определенности используем Р1С-модель частиц с сеточным ядром, которая имеет вид
1 I -л | ^^р |
Д(х - хр)=< к V1 к |х Хр|- к (17)
0, |х — хр| > к.
Рассмотрим отдельную частицу с зарядом д, которая в моменты времени ¿т = тт, ¿т+1 = (т + 1)т имеет координаты хт и хт+1.
Возможны следующие варианты перемещения частицы за временной шаг т.
А. Хг_1/2 < хт, хт+1 < хг+1/2, т.е. частица перемещается между центрами двух соседних ячеек. Тогда в соответствии с выбранной моделью (17) и формулой (15) частица в обоих положениях хт и хт+1 дает вклад только в узлы хг—1/2 и хг+1/2:
Рт—1/2 = ^(х*-1/2 - хт) = к (1 - * ^_1/2 ) , (18)
РГ+1/2 = - хт) = к (1 - Хг+1/2к- ^) . (19)
Выражения для вкладов в момент времени ¿т+1 получаются из (18), (19) простой заменой временного индекса. При этом изменение плотности в узлах сетки, связанное с перемещением частицы, выражается так:
рт_+/2 - рт—1/2 = ?^й(хг_1/2 - хт+!) - - хт) = - к^ (хт+1 - хт), (20)
РЙ^ - рт+1/2 = ^¿+1/2 - хт+Ч - дЯ(хт/2 - хт) = ^ (хт+1 - хт) . (21)
Заметим, что суммарное приращение плотности равно нулю и закон сохранения заряда выполнен.
Для выполнения схемы (14) необходимо, чтобы потоки заряда, приводящие к приращениям плотности, удовлетворяли равенствам
рт+1 т = т Лт+1/2 Зт+1/2^ (22)
Рг—1/2 - рг—1/2 = - ^ З - ^¿-1 ^ , (22)
рт+1 рт = ^ Лт+1/2 •т+1/Л (23)
Р¿+1 /2 - Р¿+1/2 = - к ^¿+1 - 3 ^ • (23)
Так как во все узлы сетки, кроме рассматриваемых двух, частица дает нулевые вклады, то при отсутствии потоков извне, последовательно используя разностное уравнение (14) во всех узлах, можно получить, что Згт+1/2 = з'т+1/2 = 0. Тогда из формул (22), (23) следуют выражения
Рт+1 _ Рт Рт+1 _ т
• т+1/2 _ 7 Р¿—1/2 Рг—1/2 -т+1/2 _ 7Рг+1/2 Рг+1/2
Зг = -к-, 3 = к-.
Подставляя в них соответствующие приращения плотности (20), (21), получаем, что из обеих формул поток определяется однозначно:
• т+1/2 _ < Жт+1 — _ ?^т+1/2
Л к т к '
Таким образом, разностное уравнение неразрывности (14) оказывается выполненным, а частица с Р1С-ядром (17) дает вклад в поток только в одном граничном узле жг, разделяющем две ячейки, между центрами которых локализована частица.
Теперь вычислим потоки заряда по формуле (16). Для определенности положим, что на всем интервале времени (¿т,£т+1) частица имеет скорость „т+1/2. Тогда ее положение в момент времени ¿т+1/2 выразится как
т жт+1 I жт
жт+1/2 _ жт + 'т+1/2 _ ж + ж
2 2 '
Пусть жг > жт+1/2. При этом в соответствии с (17) частица дает вклад в потоки в двух граничных узлах жг, жг-1, что уже противоречит выполнению разностного закона сохранения. Для потоков получаются выражения
„т+1/2 „т+1/2
Л-^1/2 = (ж - Хт+1/^ , ^-т+1/2 = (жт+1/2 - Хг-0 .
Последние, очевидно, не удовлетворяют разностной схеме (14), в которой используются приращения (22), (23).
Рассмотрим другой вариант:
Б. При жг-1/2 < жт < жг+1/2 < жт+1 < жг+3/2 частица за шаг т переходит через сеточный узел жг+1/2 слева направо. В момент ¿т частица, как и в предыдущем случае, дает вклады в узлы жг-1/2, жг+1/2 (18), (19), а в момент ¿т+1 — в узлы жг+1/2, ж^+з/2:
рт+1 _ < Л жт+1 - жг+1/Л рт+1 _ < Л жг+3/2 - ^ \ (24)
Рг+1/2 _ к ^ к ), Рг+з/2 _ к V1 к ;. (24)
В этом случае происходит изменение плотности в трех узлах сетки жг-1/2, жг+1/2, жг+3/2. Согласно разностной схеме (14), приращение плотности в узлах жг-1/2, жг+1/2 по-прежнему выражается формулами (22), (23). Приращение в узле жг+3/2 записывается аналогично:
рт+1 _ рт _ _' Лт+1/2 _ ^-т+1/Л (25)
Р г+3/2 Р г+3/2 _ к V +2 / . (25)
При этом
рт+3/2 _ рг^ _ о. (26)
Потребовав, чтобы разностный закон сохранения (14) выполнялся во всех узлах сетки,
• т+1/2 -т+1/2 п гр
в том числе и в тех, где плотность не меняется, получаем _ _ 0. Тогда из
уравнения (22) с учетом (18), (26) запишем выражение для потока в узле жг:
• т+1/2 _ к т жг+1/2 - жт
Л _ ТРг-1/2 _ к Т , (27)
а из уравнения (25), используя (24), (26),
• т+1/2 _ кЛт+1 _ < жт+1 - жг+1/2 ¡^
_ Т Рг+3/2 _ к т . (28)
Подстановка выражений (27), (28) в разностное уравнение (23) дает тождество, при этом
^^ + 1 ^^ ^ Г/ ^^ + 1 \ / ^^ \П
Р¿+1/2 - Р¿+1/2 = -^2 IIх - - (3+1/2 - х )\ •
Вычисление приращения плотности в узле ^¿+1/2 из интерполяционных формул (19), (24) приводит к тому же результату. Таким образом, для полученных формул (27), (28) разностное уравнение неразрывности оказывается выполненным. В то же время использование формулы (16) для вычисления потоков нарушает разностный закон сохранения (14).
Вариант С. Наконец, для ^¿-1/2 < хт+1 < ^¿+1/2 < хт < ^¿+з/2 при движении частицы справа налево аналогично получаем следующие выражения для плотности потоков:
.т+1/2 = к +1 = 5 х"+1 - ^¿+1/2 (29)
Л = /¿-V2 = к т ' (29)
• т+1/2 _ к _ 5 ^¿+1/2 - Хт
Л+1 = - ^¿+3/2 = к Т • ( )
Сравнивая формулы (19) и (27) - (30), можно заметить общее свойство, которое позволяет легко построить численный алгоритм. А именно, траектория частицы на шаге по времени разбивается на отдельные участки узлами ^¿-1/2, и внутри каждого участка ^¿-1/2, £¿+1/2]
т+1/2
соответствующий ток ^ подсчитывается по одной и той же формуле
• т+1/2 5 хкон хнач /01 \
* = к-Т—, (31)
где хнач и хкон — начальная и конечная координаты частицы на этом участке.
Использование полученных формул для определения плотности потока по (19), (27)-(30) в сочетании с процедурой интерполяции плотности на эйлерову сетку по формуле (15) позволяет построить консервативный алгоритм, не нарушающий разностный аналог уравнения неразрывности (13).
Заметим, что нарушение консервативности не связано с использованием конкретного способа вычисления скорости ^РТ+1/2, так как во всех формулах для вычисления потоков зарядов используются только начальные и конечные координаты частиц.
Выпишем формулы для плотности тока в двумерном и трехмерном случаях. В трехмерном случае плотности заряда и тока должны удовлетворять разностной схеме
pm+1 _ m •m+l/2 _ • rn+l/2
Pi-1/2,1-1/2,fc-1/2 pi-1/2,1-1/2,fc-1/2 _ j,i,1-1/2,fc-1/2 j,i-1,1-1/2,fc-1/2
T
• m+1/2 • m+1/2 .m+1/2 .m+1/2
+ jy,i-1/2,1,fc-1/2 — jy,i-1/2,1-1,fc-1/2 + jz,i-1/2,1-1/2,fc — jz,i-1/2,1-1/2,fc-1 (32)
где , — шаги по пространственным координатам x, y, z; , — компонен-
ты вектора плотности потока j. Дробные индексы обозначают центры ячеек, а потоки определяются в серединах граней ячеек. Как и в (14), плотности и потоки вычисляются на сдвинутых временных сетках. Таким образом, схема (32) аппроксимирует уравнение неразрывности (13) со вторым порядком 0(т2, h^). В сочетании с leapfrog scheme для лагранжевых частиц она часто используется в задачах динамики плазмы. Для выполнения уравнения (32) вклад отдельных частиц в плотность тока следует вычислять по формулам, аналогичным формуле (31), полученной для одномерного случая. Так же, как
в одномерном случае, прямолинейная траектория, соединяющая начальное (гт) и конечное (гт+1) положения частицы, разбивается на отдельные элементы, каждый из которых находится в параллелепипеде {ж*—1/2 < ж < жт/2, Уг-1/2 < у < У1+1/2, ^—1/2 < г < ^+1/2} и дает вклад в плотность тока в соответствующих этому параллелепипеду узлах сетки:
• т+1/2 _ Дж
3х,г,1—1/2,к—1/2 _ < Т
(1 - ¿„)(1 - ¿г) +
ДуДг 12кукг
. т+1/2 _ Дж
3х,г,1— 1/2,к+1/2 _ <
. т+1/2 _
3х,г,1+1/2,к—1/2 _ <
т
Дж
т
(1 - )*г -
¿у (1 - ¿г) -
ДуДг 12кукг
ДуДг 12кукг
. т+1/2 _ Дж
3ж,г,г+1/2,й+1/2 _
т
¿у ¿г +
ДуДг
у г ' 12кукг
• т+1/2 Ду
3у,г—1/2,г,£—1/2 < Т
(1 - ¿х)(1 - ¿г) +
ДжДг
12 кх кг
• т+1/2 _ Ду
3у,г—1/2,1,к+1/2 _ < Т
• т+1/2 _ Ду
3у,г+1/2,1,к—1/2 _ < Т
(1 - ¿х^г -¿х (1 - ¿г) -
ДжДг
12 кх кг
ДжДг
12к кг
• т+1/2 _ ДУ
3у,г+1/2,1,к+1/2 _
т
¿х ¿г +
ДжДг
12к кг
• т+1/2 _ Дг
3г,г—1/2,1—1/2,к _ < Т
(1 - ¿х)(1 - ¿у) +
ДжДу
12 кх ку
• т+1/2 _ Дг
3г,г— 1/2,1+1/2,к _ < Т
(1 - ¿х^у -
ДжДу
12 кх ку
• т+1/2 _ Дг
3г,г+1/2,1—1/2,к _ < Т
¿х(1 - ¿у) -
ДжДу
12 кх ку
• т+1/2 _ Дг
3у,г+1/2,1+1/2,к _ < Т
ДжДу
¿'¿у + 12кх к.
где Дж жкон жнач; Ду параллелепипеде;
_ укон - унач; Дг _ 2кон - 2нач — приращения координат в
е 1 ( жнач + жкон
¿х _ ^--ж<
х кЛ 2 '
Г 1 ( унач + укон
¿у _ кД-2--уг';
г 1 ( гнач + гкон
¿г_ к7-2--**
В двумерном случае в плоскости (х,у) формулы для вкладов в плотность тока от
частиц имеют вид
• т+1/2 _ Ах/1 X \ -т+1/2
Ах г
5-,
т
^,¿,1-1/2 = 5 т (1 ), ^,¿,1+1/2 •т+1/2 Ауп х \ -т+1/2
(1 - $х)(1 - ) +
•т+1/2 = „.т+1/2
/¿,¿-1/2,1-1/2 = У^г
т
АхАу
12 кх ку
• т+1/2 = ™,т+1/2
1/2,1+1/2 = У^г
• т+1/2 = л^т+1/2
•¿,¿+1/2,1-1/2 = У^г
(1 - -
- ) -
АхАу
12 кх ку
АхАу
12 кх ку
• т+1/2 = л^т+1/2
• ^,¿+1 /2,1+1/2 У^г
+
АхАу
12 кх ку
где ^т+1/2 — ^-компонента скорости частицы. Остальные обозначения те же, что и в трехмерном случае.
Таким образом, схема решения задачи на одном шаге разбивается на два этапа. На первом (лагранжевом) этапе по схеме (9) вычисляются скорости и координаты частиц. Здесь же определяются компоненты плотности тока j т+1/2. На втором (эйлеровом) этапе по схеме (10), (11) решаются уравнения Максвелла, т.е. определяются функции Вт+1/2, Ет+1.
3. Основные проблемы распараллеливания алгоритма
Алгоритм, реализующий метод частиц в ячейках, обладает высокой степенью внутреннего параллелизма, так как частицы непосредственно друг с другом не взаимодействуют. Несмотря на это, метод сложен для эффективной параллельной реализации. Главный вопрос, возникающий при параллельной реализации метода частиц, — каким образом распределять частицы и элементы массивов сеточных переменных между ПЭ. Следующая проблема состоит в минимизации объемов пересылок данных между ПЭ при использовании суперкомпьютеров с распределенной памятью. Например, при вычислении компонент электрических и магнитных полей необходимо знать их значения в граничных узлах сетки, которые находятся в другом процессорном элементе. Чтобы вычислить значения плотностей тока и пространственного заряда в конкретном узле, необходимо знать характеристики всех частиц, находящихся в ячейках, одной из вершин которых является данный узел, но частицы из одной или соседних ячеек могут храниться в разных процессорных элементах. Из этого следует, что пересылки избежать невозможно.
Другая проблема распараллеливания метода частиц состоит в обеспечении равномерной загрузки процессорных элементов, поскольку первоначально равномерно распределенные в пространстве частицы через несколько шагов моделирования образуют неравномерное распределение [14].
Минимальный объем пересылок во время вычисления компонент электрических и магнитных полей, плотности заряда и средней скорости достигается в случае, если характеристики частиц и значения электрических и магнитных полей в вершинах соответствующих
ячеек хранятся в памяти одного процессорного элемента, а также когда частицы, находящиеся в соседних ячейках, располагаются в одном ПЭ либо в соседних.
Универсальным способом решения этой проблемы является дублирование в памяти каждого ПЭ описания электрических и магнитных полей всего пространства моделирования. Тогда множество частиц может быть распределено между ПЭ вне зависимости от их расположения в пространстве и не требуется передавать частицы из одного ПЭ в другой в ходе моделирования. Однако при этом существенно ограничивается объем решаемой задачи, так как массивы с сеточными переменными занимают большой объем. Также требуется передавать большой объем данных между ПЭ, поскольку в каждом ПЭ находятся частицы из разных ячеек. Частицы из одной или соседних ячеек могут храниться в разных ПЭ, поэтому на каждом временном шаге необходимо суммировать значения средней плотности заряда и тока во всех узлах сетки, вычисленные в каждом ПЭ.
Предлагаемые алгоритмы распределения частиц:
1. Область моделирования разбивается поперек оси y на слои, число которых зависит от количества ПЭ (рис. 3). Таким образом, в начальный момент времени процессорные элементы имеют равномерную загрузку.
2. Каждый ПЭ хранит только необходимые ему значения сеточных переменных. Значения граничных элементов процессоры передают друг другу по кольцу.
3. На каждом временном шаге частицы, которые должны перелететь в другую ячейку, находящуюся в соседнем ПЭ, накапливаются и посылаются единым блоком. Необходимой дополнительной операцией по сравнению с однопроцессорной программой является сортировка частиц по координате для уменьшения времени пересылки.
На основании этого подхода были разработаны двумерный и трехмерный коды UMKA2D3V, UMKA3D3V, которые имеют универсальный характер, т.е. могут решать различные задачи при замене начальных и граничных условий на другие, отличные от вышеприведенных.
Результаты моделирования. Расчеты производились на CRAY T3E (операционная система UNICOS/mk 2.0.4, процессоры DEC Alpha EV5.6 с частотой 600 МГц, число процессоров — 256, Peak Performance 300 Гфлопс, две группы процессоров: 128 ПЭ по 128 Мбайт оперативной памяти и 128 ПЭ по 256 Мбайт, глобальный объем памяти — 49 Гбайт, первичный кэш 8 Кбайт, Direct Mapped, вторичный кэш 96 Кбайт, Three-way Set Associative, диск 200 Гбайт, сетевая топология 3D тор).
В качестве базового языка программирования выбран язык Fortran, компиляторы которого существуют практически на всех суперкомпьютерах. При разработке параллельной программы использовалась программная система PVM (Parallel Virtual Machine).
Результаты взаимодействия лазерного импульса с одной фольгой. Рассмотрим результаты решения задачи о взаимодействии короткого лазерного импульса с тонкой фольгой, которая расположена на расстоянии 25 длин волны лазерного импульса А
от левого края области. Амплитуда импульса а равна 5, поляризация импульса круговая. Остальные параметры таковы: длина и ширина импульса 25А и 10А соответственно. Импульс взаимодействует со слоем плазмы, толщина которого равна 0.375А, плотность соответствует а = шр/ш = 1.8. Временной шаг т составляет 0.01, число частиц в ячейке равно 8, область разделена на 320 х 320 х 320 ячеек, или 40А х 40А х 40А, всего частиц 4761600. На рис. 4 представлены плотность энергии электромагнитного поля (а), плотность электронов (б) и ионов (в) в плоскости (х, у). График распределения энергии электромагнитного поля Ш вдоль центральной прямой приведен на рис. 5, плотность электронов (сверху) и ионов (снизу) в плоскости (у, г) представлена на рис. 6.
Видно, что в результате взаимодействия создается резкий передний фронт (см. рис. 5). Плазма за время взаимодействия с лазерным лучом существенно смещается. Образующийся после взаимодействия импульс теперь можно использовать для ускорения частиц до релятивистских скоростей.
Рис. 5.
Рис. 6.
Результаты взаимодействия лазерного импульса с двумя последовательно расположенными фольгами. Рассмотрим результаты решения задачи о взаимодействии короткого лазерного импульса с двумя последовательно расположенными фольгами. Распределение полной энергии импульса при £ = 50 показано на рис. 7, а и плотности плазмы — на рис. 7, б, полной энергии импульса вдоль центральной прямой — на рис. 7, в. Первая фольга расположена на расстоянии 25А от левого края области, вторая — на расстоянии 30А. Амплитуда импульса а равна 1. Остальные параметры таковы: длина и ширина импульса равны 15А и 2.5А соответственно. Импульс взаимодействует с фольгами толщиной 0.375А, а плотность соответствует а = шр/ш = 1.1. Временной шаг т равен 0.005, число частиц в ячейке — 16.
При взаимодействии лазерного импульса с первой фольгой видно, что передний и задний фронты стали более крутыми. Часть импульса, которая прошла через первую фольгу, при взаимодействии со второй стала более крутой. Фольга при этом сместилась незначительно.
Рис. 7.
Заключение
Разработаны параллельные трехмерные и двумерные программы для численного моделирования динамики взаимодействия лазерного импульса с плазмой, в которых ячейки сетки разделены поровну между процессорными элементами, частицы имеют равномерное распределение только в начальный момент времени. В процессе счета накапливается дисбаланс частиц. Для тонкого слоя плазмы (фольги) это приемлемо, так как число частиц невелико.
В реализации программ, которые моделируют взаимодействие лазерного импульса с плазмой, заполняющей все пространство, большую часть памяти занимают характеристики частиц, а не характеристики с сеточными переменными. Нами разработан вариант, когда количество частиц в каждом ПЭ одинаково, но меняется количество сеточных слоев (код ИМКАЗБЗУа).
В дальнейшем планируется устранить недостатки первых двух вариантов и создать программу, обеспечивающую более равномерную загрузку ПЭ как по слоям сетки, так и по частицам.
Список литературы
[1] MOUROU G. A., Barty C. P. J., Perry M. D. Ultrahigh-intensity lasers: physics of the extreme on a tabletop // Phys. Today. 1998. Vol. 51. P. 22.
[2] Буллнов С. В., Вшивков В. А., ДудниковА Г. И. и др. Лазерное ускорение заряженных частиц в неоднородной плазме. Ч. I // Физика плазмы. 1997. Т. 23. С. 284-295.
[3] VSHIVKOV V. A., BULANOV S. V., NAUMOVA N. M. ET AL. Nonlinear electrodynamics of the interaction of ultra-intense laser pulses with a thin foil // Phys. of Plasmas. 1998. Vol. 5, No. 7. P. 2727-2741.
[4] Буланов С. В., Вшивков В. А., ДудниковА Г. И. и др. Взаимодействие лазерных импульсов петаваттной мощности с плазмой докритической концентрации // Физика плазмы. 1999. Т. 25, №9. C. 764-778.
[5] SENTOKU Y., LISEIKINA T. V., ESIRKEPOV T. Zh. ET AL. High density collimated beams of relativistic ions produced by petawatt laser pulses in plasmas // Phys. Review E. 200. Vol. 62, No.5. P. 7271-7281.
[6] БЕРЕЗин Ю. А., Вшивков В. А. Метод частиц в динамике разреженной плазмы. Новосибирск: Наука, 1980.
[7] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987.
[8] БэдСЕЛ Ч., ЛЕнгдон А. Физика плазмы и численное моделирование. М.: Энерго-атомиздат, 1989.
[9] LANGDON A. B., LASINSKI B. F. Electromagnetic and relativistic plasma simulation models // Meth. Comput. Phys. 1976. Vol. 16. P. 327-366.
[10] BORIS J. P. Relativistic plasma simulation - optimization of a hybrid code // Proc. Fourth Conf. Num. Sim. Plasmas. Washington, 1970. P. 3-67.
[11] MARDER B. A method for incorporating Gauss' law into electromagnetic PIC codes //J. Comput. Phys. 1987. Vol. 68. P. 48-55.
[12] VILLASENOR J., BUNEMAN O. Rigorous charge conservation for local electromagnetic field solver // Computer Phys. Comm. 1992. Vol. 69. P. 306-316.
[13] Вшивков В. А., Вшивков К. В., ДудниковА Г. И. Параллельный 3D алгоритм решения задачи взаимодействия лазерного импульса с плазмой / / XVI школа-семинар "Информ. технологии в задачах мат. моделирования": Тез. докл. (http://www.ict.nsc.ru/ws/ct-2000/).
[14] Вшивков В. А., Краева М. А., МАЛЫшкин В. Э. О реализации метода частиц на мультипроцессорах. Новосибирск, 1995 (Препр. СО РАН ВЦ; №1052).
Поступила в редакцию 4 января 2001 г., в переработанном виде —15 января 2001 г.
