CC BY
17
УДК 533.9
Ю.А. Кропотина, М.Ю. Густов, A.M. Красильщиков, К.П. Левенфиш, Г.Г. Павлов
МНОГОПРОЦЕССОРНЫЙ ТРЕХМЕРНЫЙ ГИБРИДНЫЙ КОД ДЛЯ МОДЕЛИРОВАНИЯ МИКРОСКОПИЧЕСКИХ ЯВЛЕНИЙ В КОСМИЧЕСКОЙ ПЛАЗМЕ
Теоретическое исследование бесстолкнови-тельных процессов в космической плазме представляет большой интерес для современной астрофизики, поскольку такие процессы играют большую роль в динамике источников с экстремальным энерговыделением, в частности, остатков сверхновых звезд, пульсарныхтуманностей и активных галактических ядер. Поскольку задачи, связанные с изучением микроскопических процессов в плазме, являются существенно нелинейными, аналитические подходы не позволяют далеко продвинуться в их решении. С другой стороны, численное моделирование, основанное на фундаментальных физических законах, позволяет давать как качественные , так и количественные интерпретации современных многоволновых наблюдений космических объектов с экстремальным энерговыделением, и оценивать параметры этих объектов.
Конкретная методика численного моделирования выбирается в зависимости от пространственных и временных масштабов исследуемых процессов [1]. В случае, когда динамика плазменной системы в основном определяется ионами, пространственные и временные масштабы для которых существенно больше, чем для электронов, одним из эффективных методов является гибридное моделирование, в котором ионы описываются как частицы, а электроны представлены как безмассовая нейтрализующая жидкость. Такой подход расширяет область параметров, в которой можно производить моделирование, и незаменим на данном этапе развития вычислительной техники для исследования крупномасштабных процессов. Следует отметить, что экономия вычислительных ресурсов за счет использования пространственно одномерных или двумерных кодов с учетом электронов в виде частиц ведет к возникновению нефизических эффектов [2,3] и не может быть признана адекватной.
Однако даже в приближении жидких электронов вычислительные ресурсы, требуемые для создания трехмерных моделей, превышают возможности современных персональных компьютеров, в связи с чем необходимо создание кодов, способных работать на многопроцессорных кластерах.
В данной работе представлен трехмерный гибридный код второго порядка точности по времени и координате, адаптированный для работы на многопроцессорных вычислительных системах, и дан пример его использования для моделирования резонансной ион-ионной неустойчивости в космической плазме.
Система уравнений для микроскопического описания бесстолкновительной плазмы
В целях уменьшения погрешности, связанной с машинной точностью, рационально использовать безразмерные физические величины, для чего выбрана следующая нормировка [4, 5]:
? = г-П; (1)
г = гЩ\ (2)
в = В/В0; (3)
Р = Р/Ро, (4)
где В0, р0 — индукция магнитного поля и плотность в невозмущенной плазме (в случае моделирования ударной волны нормировочные величины берутся перед фронтом), Х = еВ0/тс циклотронная частота, /,. =с/ю — инерционная длина иона, ю = (4лр0е2/т2)1^2 — плазменная частота, е, т — заряд и масса протона, с — скорость света.
В подавляющем большинстве астрофизических задач исследуется плазма, преимущественно состоящая из водорода, поэтому здесь приведены нормировочные параметры для водородной
плазмы. Из нормировок (1) и (2) следует, что скорости измеряются в альвеновских скоростях, т. е. вместо скорости налетающего потока задается число Маха:
V = У1Уа=МаУ
где =
Заряды и массы частиц измеряются в зарядах и массах протона. Однако представленный код позволяет моделировать также и многокомпонентную плазму, при этом для нормировки используются взвешенные величины, а вид уравнений не меняется.
Поскольку ионы описываются как частицы, для них используется уравнение Лоренца движения в электромагнитном поле. К этим уравнениям добавляются уравнения Максвелла для вычисления полей и уравнение Эйлера для описания движения электронов. В приведенной нормировке система уравнений, определяющая поведение системы, может быть записана в следующем виде [6]:
dr,
dt Z i
dt тк
VxB = J;
dVk
— = -VxE;
dt
' dt
= -«„E + J„xB.
(5)
(6)
(7)
(8)
(9)
J = i¡ +Je;
J;
к
Jк -{na)k,
(10) (11)
(12)
где I,., Зк — токи ионов и частиц сорта к\ (пи)к — среднее значение произведения скорости на концентрацию для данного сорта частиц.
Если пренебречь массой электрона по сравнению с массой протона и учесть квазинейтральность
то уравнение (9) можно переписать в следующем виде:
Е = -(УхВ)хВ--(1/хВ). (13)
п п
Для дальнейших вычислений удобнее привести его к виду
1 dPfk
Ек= —(J,xB)fc---
1 —i
п
п dXj
(14)
В системе (5)—(9) гь \к — координаты и скорости А;-го сорта частиц. 2к, тк — заряд и масса в единицах заряда и массы протона, Е, В — электрическое и магнитное поля, и^, пе — скорость потока и концентрация электронов соответственно; Зе — ток электронов, ,1 — электрический ток. Они вычисляются по следующим формулам:
В
P?k = —bjk - BjBk ления.
Алгоритм работы гибридного кода
Каквидно из уравнений (5) — (9), поля ича-стицы взаимно зависимы, что делает задачу самосогласованной. Задачи такого рода удобно решать с помощью известного метода липфрог («leapfrog»), имеющего второй порядок точности по времени. При этом значения электромагнитного поля известны на четных шагах, а положения частиц — на нечетных.
Структура уравнений для вычисления электромагнитных полей аналогична структуре уравнений магнитной гидродинамики (МГД), вследствие чего используются сходные методы решения. Для их реализации вводится декартова пространственная сетка. На гранях ячеек задаются нормальные значения магнитного поля, для вычисления приращения которых используются значения электрического поля на ребрах сетки [7].
Основная последовательность операций после инициализации начальных параметров соответствует стандартному алгоритму работы гибридного кода [4—6]. Циклически поочередно выполняются следующие процедуры.
1. Вычисление моментов (Moment collector). На основе текущих положений и скоростей частиц производится подсчет плотностей и токов во всех ячейках по формулам (11), (12).
2. Вычисление электромагнитных полей (Field solver). С использованием уравнений (7), (8), (14) вычисляются электрические и магнитные поля в ячейках.
3. Перемещение частиц (Particle mover). Производится перемещение и изменение скоростей всех частиц согласно уравнениям (5), (6).
Код реализован на языке С++ с использованием методов объектно-ориентированного программирования.
Вычисление моментов. При вычислении моментов в каждой ячейке производится суммирование по всем находящимся в ней частицам. Эффект дискретизации приводит к возникновению численных шумов на этом шаге. Распространенным методом их подавления является введение макрочастиц — частиц конечного размера [8, 9]. Однако тесты показали, что применение метода макрочастиц практически не меняет уровень флуктуаций. При этом ресурсоемкость увеличивается в восемь раз (для трехмерного кода). В связи с этим в коде используется точечное представление частиц.
Алгоритм вычисления электромагнитных полей. В коде используется алгоритм вычисления магнитного поля с сохранением нулевой дивергенции, основанный на методе Годунова для решения уравнений МГД[7, 10, 11]. Он включает параболическую реконструкцию магнитного поля внутри ячейки [7] и линеаризованную задачу о распаде разрыва (Riemann solver) [10,12].
Для реализации трехмерного бездивергентного алгоритма вводятся значения нормальной компоненты магнитного поля на гранях и значения электрического поля на ребрах. Считается, что центр ячейки имеет целые координаты, а грани — полуцелые.
Для повышения точности при вычислении электрического поля, а также приращений магнитного используется дополнительное разбиение по времени. За время одного шага по времени At с неизменным положением частиц и неизменной конфигурацией токов, п раз выполняется алго-
dt = 9/n
ВЦ=В(0;
В{=В'0-Л-(УхЕ'0), Bj=5¿-fiMVxEj), В (t + dt ) = 0,5(В 0 + В2),
В(г + Д0 = В( t + ndt).
На каждом шаге dt решается задача Римана для расчета электрического поля на ребрах, а затем вычисляется магнитное поле на гранях. При этом выполняются следующие действия.
1. С использованием известных значений нормальной компоненты магнитного поля на гранях строится реконструкция поля внутри ячейки [7]:
Вх(х,у,г) =
= а0 + ахх + ауу + агг + ашх2 + а^ху + ахгхг;
Ву(Х,у,2) =
= Ь0 +Ьхх + Ьуу + Ьгг + Ьхуху + Ьууу2 + Ьугуг,
=
= со +Схх + суу + сгг + схгхг + сугуг + сггг2.( 13)
Для выполнения условия сохранения нулевой дивергенции магнитного поля коэффициенты разложения рассчитываются следующим образом:
./+1/2 _ Ы-1/2)ч
(Bl^/¿-Blx:i/¿>)
9х
av = — у 2
а, =-
9УВХ
,/+1/2
9УВХ
-1/2Л
AJ
Aj
А Л
i+1/2
А
/-1/2'
А
А
йху=~Ах
Аув?1/2 АуВ^/2
АУ
АУ
А
АВМ/2 А В;'-1/2Л
1—У 7 ХУ .. ¿—Ж ..
А
1
А
ахх 2 —(Ьху +СхгУ,
% =
(В^2+В^2)
-аг
9х
Система обозначений пояснена на рис. 1, изображающем ячейку с координатами (/,у, к).
Рис. 1. Ячейка декартовой сетки. Магнитные поля заданы на гранях, электрические — на ребрах сетки; обозначены точки, в которых решается задача Римана для нахождения значений электрического поля
Реконструкция производится с использованием лимитера тшпюс!:
= minmod^' где
minmod(a,6) =
9у
'+1/2,7+1 д'+1/2,j _ R/+1A7_1
ж
-к
fsign(a)min(|a|,|6|), еслиа6>0;
[О,
если ab < 0.
шаг 1 выполняется так же, как и на промежуточной стадии;
на шаге 2 задача Римана решается не для ребер, а для граней ячеек; при этом находятся компоненты тензора магнитного давления P¡j. Для достижения второго порядка точности по координате используется полиномиальное разложение Ру на грани, и задача Римана решается в четырех точках;
на шаге 3 электрическое поле в ячейке вычисляется путем численного дифференцирования значений Ру на гранях. Таким образом, на шаге Field solver считается только градиентная часть электрического поля (второе слагаемое в правой части уравнения (14)). Часть, связанная с током ионов, рассмотрена в следующем разделе.
Магнитное поле в ячейке считается путем усреднения реконструкции (13). При этом среднее значение компоненты поля следует выражению [7]:
_ (в•
Bf _-
M/2+Bi-i/2
l)
2. Для нахождения значения электрического поля на ребрах решается задача Римана о распаде разрыва, возникающего на границах ячеек после реконструкции поля. Для достижения второго порядка точности значение электрического поля на ребре аппроксимируется полиномом третьей степени по координате. Это приближение позволяет представить среднее электрическое поле каждого ребра как полусумму значений поля в двух его точках, обозначенных на рис. 1. Поскольку ребро является совмещенной границей четырех ячеек, задача о распаде разрыва в обеих точках решается дважды (на каждой прилегающей грани), а результат усредняется.
3. С использованием полученного значения электрического поля на ребрах, из уравнения (8) находится приращение магнитного поля на гранях.
По окончании цикла «leapfrog» выполняется финальная стадия алгоритма, отличающаяся от промежуточных шагов:
_
- 2 "йхх 6 '
Движение частиц. Как отмечено выше, в коде используется метод «leapfrog», т. е. для осуществления перемещения частиц t ^ t + At используются значения электрического и магнитного поля в момент времени í + 0,591. Значения магнитного поля и градиентной части электрического поля определяются на шаге Field solver. Однако для корректного учета в уравнении (14) слагаемого, зависящего оттока ионов, необходимо аппроксимировать его изменение в течение шага по времени, вызванное движением частиц. Для этого на шаге Particle mover производится линейная экстраполяция тока ионов и в линейном приближении учитывается изменение электрического поля на протяжении шага по времени, что позволяет сохранить второй порядок точности по времени при решении самосогласованной задачи и избежать накопления ошибки, приводящего к нефизическим эффектам. Здесь следует отметить, что линейная экстраполяция динамики электрического поля в течение шага по времени некорректна в случае введения частиц конечного размера.
Перемещения частиц вычисляются при помощи разложения формул (5), (6) в ряд Тейлора, что упрощает учет динамики электрического поля. Разложение осуществляется до четвертого порядка во избежание внесения в вычисле-
ния дополнительной погрешности на шаге Particle mover.
Применимость использованного метода. Метод «leapfrog», реконструкция магнитного поля, полиномиальная аппроксимация электрического поля при решении задачи Римана и самосогласованный учет динамики электрического поля при движении частиц обеспечивают второй порядок точности алгоритма по времени и координате.
Для проверки адекватности результатов моделирования необходимо контролировать выполнение фундаментальных законов сохранения в течение работы кода. Как было показано выше, алгоритм вычисления полей обеспечивает точное сохранение нулевой дивергенции магнитного поля.
Тесты показали выполнение законов сохранения энергии и импульса с точностью второго порядка по времени при корректно выбранных параметрах дискретизации, как показано ниже.
Как отмечено выше, в коде используется декартова система координат с равномерным разбиением. Размер ячейки, в зависимости от поставленной задачи, может варьироваться от 0,1// до нескольких инерционных длин иона [4]. Нижний предел связан с корректностью приближения безмассовых электронов, поскольку на мелких масштабах становятся существенными электронные эффекты. Верхний предел определяется минимальными характерными длинами исследуемых процессов. Шаг по времени не должен превышать 0,2 для хорошего разрешения траектории движения иона в магнитном поле [4]. Кроме того, соотношение временного и пространственного шагов должно быть таким, чтобы частица за один шаг по времени проходила расстояние меньше размера ячейки (критерий Куранта). Во многих задачах критерий Куранта оказывается жестче, и шаг по времени определяется им.
Ограниченное число частиц в ячейке приводит к появлению численного шума при вычислении электрических и магнитных полей. Численный шум возникает не только при исследовании неравновесных процессов, но и в простом тестовом примере — моделировании однородно распределенного водорода в постоянном магнитном поле. При этом наблюдается генерация флуктуаций поля, уровень которых хорошо описывается зависимостью
где N — число частиц в ячейке.
Такая зависимость указывает на статистическую природу флуктуаций, связанную со случайным распределением частиц в начальный момент времени. Эффективный способ подавления этих флуктуаций — увеличение числа частиц в ячейке, однако это ведет к уменьшению статистического веса частиц и увеличению погрешности, связанной с точностью машинного числа, а также к существенному росту потребляемых машинных ресурсов.
В большинстве решаемых задач число частиц в ячейке составляло 100—200, что обеспечивало среднеквадратичный уровень флуктуаций не более 10 %. Фактическое число частиц в физических задачах на много порядков больше, чем число частиц, которое возможно реализовать в программе, поэтому каждая частица в ячейке может быть интерпретирована как совокупность множества ионов и при вычислении моментов она суммируется с определенном весом.
Распараллеливание вычислений. В связи с большой ресурсоемкостью трехмерного алгоритма код адаптирован для многопроцессорного моделирования с использованием технологии MPI (обмен сообщениями между процессорами) [13]. Поскольку на практике часто приходится моделировать области с размером вдоль одной координаты, существенно превышающим размеры по двум другим, распараллеливание программы было осуществлено путем разделения пространства моделирования между несколькими процессорами и введения нового типа граничных условий — границы между областями, обрабатываемыми в различных процессорах. На этой границе на каждом шаге вычислений происходит обмен значениями величин поля и потоков частиц между процессорами, обрабатывающими две соседние области. Для реализации процедуры обмена создаются вспомогательные ячейки пространственной сетки.
Использование технологий параллельного программирования позволяет на порядки увеличить физический размер области моделирования, точность и скорость расчета. Следует отметить, что разделение области вычислений между процессорами рационально только в случае, когда длина, приходящаяся на один процессор, составляет как минимум 3 ячейки. В противном случае ресурсы памяти, необходимые для реализации
вспомогательных ячеек, превысят ресурсы, потребляемые вычислительной частью.
Пример использования кода для моделирования
резонансной ион-ионной неустойчивости
В работе [4] приведен пример моделирования задачи о возникновении резонансной электромагнитной ион-ионной неустойчивости. Ниже приведены результаты аналогичных расчетов для демонстрации корректной работы кода на примере аналитически решаемой задачи. Резонансная ион-ионная неустойчивость возникает при прохождении слабого (я, << я0) пучка ионов через равновесную фоновую плазму п0 =п1+пх) вдоль однородного магнитного поля В()х с относительной скоростью У{ >> ¥а. Взаимодействие протонного пучка с протонной фоновой плазмой рассматривается в системе покоя нейтрализующих электронов. В ней пучок движется со скоро-
стью У{ (1 - я, /я0), а фоновая плазма — со скоростью ^ =- (^ п0)У|. В таких условиях возникают нарастающие резонансные циклотронные колебания. В работе [4] при помощи одномерного гибридного кода смоделирована неустойчивость плазмы со следующими параметрами: скорость пучка относительно фоновой ==
ношение теплового к магнитному давлению Р = /В0 = 1 для пучка и для фоновой плазмы. Трехмерные расчеты с теми же параметрами привели к результатам, согласующимся с приведенными в работе [4].
На рис. 2 показаны результаты моделирования — фазовое пространство, а также профили магнитного поля и плотности пучка. Все величины приведены в нормировке по формулам (1)—(4). Вывод данных был осуществлен в момент времени г = 40П-1. К этому моменту произошло суще-
К
5
20
40
60
80
Рис. 2. Результат трехмерного гибридного моделирования резонансной неустойчивости: а — фазовое пространство х— Уу пучка; б, в — профили магнитного поля и плотности пучка в момент времени / =
ственное нарастание гармоник с длиной волны, близкой к 50/;, что согласуется с теоретическими предсказаниями. Соответствующие гармоники наблюдаются как для скорости, таки для профилей плотности пучка и поля, причем отчетливо видна разность фаз между флуктуациями поперечной скорости пучка и флуктуациями профилей плотности и поля, составляющая около 90°. На больших временах неустойчивость приводит ктермализации пучка.
Итак, разработан трехмерный многопроцессорный гибридный код для описания астрофизической плазмы с точным сохранением нулевой дивергенции магнитного поля и выполнением законов сохранения энергии и импульса с точностью второго порядка по времени. Разработанный алгоритм включает метод Годунова с реконструкцией магнитного поля и решением задачи Римана, традиционно используемый в магнитогидродинамических кодах и впервые адаптированный для гибридного кода. В алго-
СПИСОК J
1. Winske, D. A nonspecialist's guide to kinetic simulations of space plasmas [Text] / D. Winske, N. Omidi // Journal of geophysical research.— Vol. 101,- А8,— 1996,- P. 12287-17303.
2. Giaealone, J. Charged-particle motion in multidimensional magnetic-field turbulence [Text] /J. Giaealone, J.R. Jokipii // Astrophys. Journ. 1994.— Mo 430,- P. 137-140.
3. Jones, F.C. Charged-particle motion in electromagnetic fields having at least one ignorable spatial coordinate |Text| / F.C. Jones, J.R. Jokipii, M.G. Baring // Astrophys. Journ.- 1998,- № 509,- P. 238243.
4. Winske, D. Hybrid codes. Methods and applications. Computer space plasma physics: Simulation techniques and softwares [Text] / D. Winske, N. Omidi, edited by H. Matsumoto, Y. Omura // Tokyo: Terra Scientific, 1993,- P. 103-160.
5. Filippyehev, D.S. Hybrid simulation of space plasmas: models with massless fluid representation of electrons. 1. Collisionless shocks [Text] / D.S. Filippyehev // Computational Mathematics and Modeling.- 2000,- Vol. 11,- № 1,- P. 15-39.
6. Matthews, A.P. Current advance method and cyclic leapfrog for 2D multispecies hybrid plasma simulations I Text I / A. P. Matthews// Journ. Сотр. Phys.- 1994,- № 112. P. 102-116.
7. Balsara, D.S. Second-order-accurate schemes
ритме не применяется традиционно используемый метод частиц конечного размера, что позволяет самосогласованно вычислять динамику электрического поля при перемещении частиц. Построенный код протестирован путем сравнения с одномерным аналогом и аналитическими результатами для упрощенного случая резонансной ион-ионной неустойчивости и может быть использован для моделирования бес-столкновительных процессов в астрофизической плазме.
Авторы выражают благодарность профессору A.M. Быкову за постановку задачи и неоценимую помощь в ее выполнении, а также Московскому и Санкт-Петербургскому отделениям МСЦ РАН за предоставление вычислительных ресурсов (кластер МВС-ЮОК МСЦ РАН и кластер HP).
Работа поддержана Министерством образования и науки Российской Федерации (договор 11.G34.31.0001 с СПбГПУ и ведущим ученым Г.Г. Павловым), а также грантом РФФИ 11-02-00429 и Фондом Дмитрия Зимина «Династия».
for magnetohydrodynamics with divergence-free reconstruction |Text| / D.S. Balsara // Astrophys. Journ. Suppl, Sen- 2004,- № 151,- P. 149-184.
8. Pritchett, P.L. Particle-in-cell simulation of plasmas— a tutorial I Text I / P.L. Pritchett // Space plasma simulation: Proc. of the Sixth International School/ Symposium, 1SSS-6, Garching, Germany, 3-7 September, 2001. Edited by J. Büchner, С.Т. Dum, M. Scholen Berlin: Schaltungsdienst Lange o.H.G.— 2001,- P. 1.
9. Dawson, J.M. Particle simulation of plasmas I Text I / J.M. Dawson // Rev. Mod. Physics.- 1983.— Vol. 55,- No. 2,- P. 403-447.
10. Годунов, С.К. Численное решение многомерных задач газовой динамики [Текст]/ С.К. Годунов,— М.: Наука, 1976,— 400 с.
11. Gardiner Т.А. An unsplit Godunov method for ideal MHD via constrained transport [Электронный ресурс I / Т. A. Gardiner, J.M. Stone // arxiv:astro-ph/ 0501557 vl. 25.01.2005
12. Balsara, D.S. Linearized formulation of the Riemann problem for adiabatic and isothermal magnetohydrodynamics [Text| / D.S. Balsara // Astrophys. Journ. Suppl Ser..— 1998— N° 116.— R 119-131.
13. Воеводин, B.B. Параллельные вычисления |Текст| / B.B. Воеводин, Вл. В. Воеводин,— СПб.: БХВ-Санкт-Петербург, 2002,- 608 с.
