CC BY
8Научный вестник НГТУ. - 2010. - № 3(40)
УДК 519.63
Вычисление на суперЭВМ температуры
при моделировании высокотемпературной плазмы
методом частиц в ячейках*
В.А. ВШИВКОВ, А.В. СНЫТНИКОВ
Проводится моделирование теплопроводности в высокотемпературной плазме в трехмерной постановке. Модель построена на основе метода частиц в ячейках. Параллельная реализация выполнена с помощью метода декомпозиции области: используется смешанная эйлерово-лагранжева декомпозиция. Наиболее сложной задачей при этом является вычисление температуры плазмы. В работе предложена методика вычисления температуры при моделировании плазмы методом частиц и доказана ее корректность.
Ключевые слова: математическое моделирование, плазменные неустойчивости, метод частиц в ячейках, параллельные вычисления
ВВЕДЕНИЕ
Метод частиц в ячейках получил широкое распространение для моделирования нестационарных задач физики разреженной плазмы в связи с бурным развитием вычислительной техники, в том числе появлением многопроцессорных комплексов. К достоинствам метода относится внутренний параллелизм: траектории модельных частиц могут вычисляться независимо друг от друга, в том числе и на разных процессорах [1]. Фактически это единственный метод для моделирования упомянутых выше задач [2].
В методе частиц каждая частица становится носителем некоторого набора характеристик среды, таких как заряд, масса, импульс, кинетическая энергия и т. д. Для того чтобы связать с каждой модельной частицей определенную температуру, необходимо вычисление температуры по ансамблю частиц, но при этом возникает проблема отделения температуры от нефизических эффектов (шумов), возникающих в методе частиц. Источниками нефизических шумов в методе частиц в ячейках являются дискретность модели - наличие сетки, на которой вычисляются распределения плотности, скорости, тока, малое число модельных частиц, большой временной шаг. Обычно при решении задач физики плазмы с помощью метода частиц в ячейках температура плазмы считалась равной нулю, т. е. разброс частиц по скоростям отсутствовал. С течением времени в плазме развивались счетные шумы и расчеты можно было считать верными до тех пор, пока влияние счетных шумов было незначительным. Единого подхода к решению проблемы нефизических шумов не выработано, нет даже общепринятого количественного определения, что такое нефизический шум в методе частиц. Для уменьшения уровня нефизических шумов чаще всего либо увеличивают количество частиц, что не всегда возможно в силу ограниченности ресурсов ЭВМ, либо модифицируют форму частицы. Трудность заключается в том, что при этом меняется и значение модельной температуры.
С развитием вычислительной техники появилась возможность решать более сложные задачи, моделировать более тонкие физические эффекты. Но при решении задач с большим чис-
* Получена 24 мая 2010 г.
Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований, гранты 08-01-615 и 08-01-622, а также интеграционных проектов СО РАН № 103, № 113 и № 26, и при поддержке гранта Президента РФ МК-3562.2009.9
лом частиц и в течение большого времени стали проявляться недостатки метода, которые раньше не играли существенной роли. Эти недостатки связаны, прежде всего, с проявлением счетных шумов.
Актуальность настоящей работы связана с тем, что в экспериментах на многопробочной магнитной ловушке ГОЛ-3 (ИЯФ СО РАН) наблюдается понижение электронной теплопроводности на 2-3 порядка по сравнению с классическим значением [3]. При этом происходит нагрев плазмы до термоядерных параметров и появляется необходимость моделирования температуры.
Плазма в экспериментах на ГОЛ-3 является существенно неравновесной, что делает невозможным применение численных методов гидродинамического типа для моделирования этой плазмы. Поэтому для решения задачи об аномальной теплопроводности в данной работе используется метод частиц в ячейках. Основная сложность заключается в том, что в методе частиц в ячейках не существует определения температуры, вдобавок требуется уточнить, что такое температура при отсутствии гауссова распределения частиц по скоростям (т. е. в неравновесной плазме).
1. ОПИСАНИЕ МОДЕЛИ
Численная модель высокотемпературной бесстолкновительной плазмы включает уравнение Власова и систему уравнений Максвелла:
1Г+У1Г+^ = = ** [Е+-[у^;
„ 4я . 1 Ж rot B = — j +--;
c c dt
1 dB
rotE =---; .
c dt
divE = 4лр; divB = 0.
В данной работе используется процедура решения этих уравнений, описанная в [4]. Далее все уравнения буду приводится в безразмерном виде. Для обезразмеривания используются следующие базовые величины: скорость света c = 3х1010 см/с; плотность плазмы n0 = 1014 см3; плазменная электронная частота юе = 5.6*104 n12 с-1.
Уравнение Власова решается методом частиц в ячейках. В этом методе вся плазма моделируется набором отдельных частиц, каждая из которых характеризует движение многих физических частиц. Характеристики уравнения Власова описывают траектории этих частиц. Уравнения этих характеристик имеют вид
f = -( E + [ v. B]),
dr _1
— = V. p = yv. у dt
_1=Vi_
2
Здесь коэффициент к для электронов равен 1, а для ионов - отношению масс электрона и иона.
Для нахождения электрических и магнитных полей используется схема, в которой поля определяются из разностных аналогов законов Фарадея и Ампера [4]:
т.т+1/2 т.т-1/2
B-—-=rot hEm,
X
m +1 m
E - E = -4*jm+1/2 + rot hBm+1/2.
X
Эта схема имеет второй порядок аппроксимации по пространству и по времени.
Задача ставится следующим образом. В начальный момент в трехмерной области решения (размер области L), которая имеет форму прямоугольного параллелепипеда
0 < х < L, 0 < y < L, 0 < z < L,
находится плазма, состоящая из электронов и ионов. Частицы (как реальные, так и соответствующие им модельные) распределены по области равномерно. Задается плотность плазмы и температура электронов, температура ионов считается нулевой. Масса иона равна 1836 масс электрона (ион водорода), заряд иона равен +1. Дополнительно в области присутствуют электроны пучка, которые также распределены по области равномерно (предполагается, что пучок уже вошел в моделируемую область). Электроны пучка отличаются от электронов плазмы тем, что не имеют разброса по скоростям (их температура равна 0) и их энергия на несколько порядков выше. Модельные электроны пучка имеют меньшую массу, нежели модельные электроны плазмы (отношение их масс равно отношению плотности плазмы и плотности пучка).
Итак, исходными параметрами задачи являются плотность n и температура T плазмы, отношение v плотности плазмы и плотности пучка, энергия пучка и масса иона. Требуется найти распределение ионов и электронов по энергиям и инкременты плазменных неустойчивостей.
2. ВЫЧИСЛЕНИЕ ТЕМПЕРАТУРЫ ПО РАСПРЕДЕЛЕНИЮ МОДЕЛЬНЫХ ЧАСТИЦ В МЕТОДЕ ЧАСТИЦ В ЯЧЕЙКАХ
Для численного моделирования неравновесной плазмы (т. е. плазмы при наличии каких-либо неустойчивостей) метод частиц в ячейках является наилучшим. В то же время в методе частиц в ячейках существует проблема нефизических эффектов - так называемых шумов. Особенно это актуально при вычислении температуры, так как приходится определять, является ли разброс по скоростям модельных частиц физическим (температура) или счетным (шумы) эффектом.
Основной источник нефизических шумов в методе частиц в ячейках - наличие сетки, на которой вычисляются распределения плотности, скорости, тока. Поэтому при расчете температуры в методе частиц в ячейках необходимо определить, насколько вычисляемая по частицам температура подвержена влиянию сетки и как изменяются значения температуры при увеличении числа узлов сетки при уменьшении пространственного шага. В настоящей работе были получены ответы на эти вопросы.
В настоящей работе предложена методика для вычисления температуры по распределению модельных частиц в методе частиц в ячейках. Температура рассчитывается как дисперсия импульсов модельных частиц:
1 Н 2 1 Н Т = Ж § & - Р0 } , Р0 = N § р .
При этом в вычислении температуры не используются массы частиц, что было бы не вполне корректным для релятивистских частиц. Кроме того, температура оказывается средней кинетической энергией частиц в соответствии с физическим определением [5].
На рис. 1 показана средняя (линия с маркерами) и минимальная по ячейкам сетки температуры электронов.
Средняя температура Минимальная температура
100 150
Число частиц в ячейке
Рис. 1. Температура в зависимости от числа модельных частиц
Для того чтобы вычисляемое таким образом значение температуры могло считаться физическим, необходимо показать, что оно не зависит от вычислительных параметров, в первую очередь числа модельных частиц в ячейке N и шага сетки И (или, что то же самое, числа узлов сетки М при фиксированной длине области L, И = L/M).
Для расчетов на трехмерной сетке было показано, что вычисляемая по частицам температура не зависит от числа модельных частиц в ячейке Ыс, если это число достаточно велико. Из рис. 1 видно, что температура незначительно изменяется при изменении Ыс.
Зависимость температуры от числа узлов сетки исследовалась на примере задачи о выравнивании электронной и ионной температур в плазме. Изначально электроны имели некоторую температуру (соответствующую 56 эВ), для ионов была задана нулевая температура. Так как начальные значения модельной температуры на разных сетках при всех тех же параметрах немного различаются, для сравнения использовалось относительное изменение температуры которое определяется как отношение изменения температуры электронов дТ к начальному значению температуры электронов Т0: = 8Т/Т0.
Как видно из рис. 2, величина зависит от шага сетки (числа узлов). Тем не менее это сходящая зависимость.
Выводы относительно зависимости температуры от шага сетки и числа частиц находятся в согласии с теоретическими оценками [2, 6].
Исследовался также нефизический разогрев плазмы вследствие флуктуаций самосогласованного электрического поля. На рис. 3 показана зависимость от шага сетки И отношения сеточной температуры электронов после завершения расчета (400 тыс. временных шагов, 100 частиц в ячейке) к начальному значению сеточной температуры электронов. Отдельно показаны максимальные и минимальные (по всем ячейкам сетки) значения температуры. Показанные на рисунке значения шага сетки соответствуют числу узлов 8, 16, 32 и 64. Начальное значение температуры (дисперсии импульсов) 0.0001.
На рис. 3 приведено максимальное значение отношения конечной температуры к начальной. Это означает, что лишь в отдельных ячейках сетки температура возрастает в 15-20 раз. \
бе-05 -
5е-05
1е-05 -
Число узлов сетки
Рис. 2. Изменение температуры электронов в зависимости от шага сетки
h
Рис. 3. Саморазогрев модельной плазмы. Отношение максимальной и минимальной сеточной температуры после завершения расчета к начальному значению
Саморазогрев плазмы под влиянием флуктуаций - это, очевидно, нефизическое явление. То что саморазогрев не зависит от шага сетки, никак не характеризует правильность расчета в целом. В данном случае ставится цель показать, что предложенная методика вычисления не зависит от шага сетки, как это и показано на рис. 3.
Наконец, показано, что коротковолновые гармоники температуры в трехмерном случае содержат лишь незначительную часть полной энергии спектра. Вычислялось отношение суммы квадратов коротковолновых гармоник к сумме квадратов всех гармоник спектра:
R У } = {п} I I I р ) .
i=0.9MX к=0.9МГ I=0.9MZ
Здесь Мх, Му, М2 - число узлов сетки по координатам X, У и 2 соответственно, Р(П)цк - фурье-гармоника сеточной величины п с волновыми номерами /, I, к в момент времени t а S{n} определяется следующим выражением:
М X Му М 2
S {п }=1 II 1р (П,к,/ )2.
г=1 к=1 I=1
Вычисленная таким образом величина R не превышала 0.1 % как для температуры (R{T}), так и для плотности электронов (R{n}):
0.0015 < R{T} < 0.004,
0.00062 < R{n} < 0.00078.
Именно в этих величинах вследствие небольшой массы электронов наиболее ярко проявляются нефизические эффекты. Если бы коротковолновые гармоники преобладали, т. е. содержали бы более 10 % энергии спектра (R{n} > 0.1), то это означало бы, что большая часть взаимодействий происходит на небольших масштабах (1-2 шага сетки), что стохастические перемещения модельных частиц внутри ячейки вносят существенный вклад в динамику частиц и в целом результат вычисления температуры определяется сеткой.
Результаты проведенных расчетов свидетельствуют, что для достаточно подробных сеток температура, определяемая по распределению модельных частиц, не зависит от шага сетки и числа модельных частиц. Кроме того, показано, что сетка не оказывает определяющего влияния на результат расчета, так как относительная плотность коротковолновых сеточных гармоник имеет незначительную величину. Можно еще добавить, что величина отклонения полной энергии системы от начального значения не превышает 1 %, а полного импульса - 0.1 % за все время расчета. Таким образом, можно утверждать, что температура, вычисленная как дисперсия импульсов модельных, частиц имеет физический смысл.
3. ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ
Необходимость параллельной реализации диктуется потребностью проводить расчеты с большим количеством модельных частиц и на больших сетках. Параметры типичного расчета: 2 млн узлов сетки, 300 млн модельных частиц. Распараллеливание выполнено методом декомпозиции расчетной области по направлению, перпендикулярному направлению движения электронного пучка. Используется смешанная эйлерово-лагранжева декомпозиция.
Сетка, на которой решаются уравнения Максвелла, разделена на одинаковые подобласти по одной из координат (^ подобластей). С каждой подобластью связана группа процессоров (в том случае, когда вычисления производятся на многоядерных процессорах, процессором для единообразия будет именоваться отдельное ядро). В каждой группе КТ процессоров, т. е. всего К*КТ процессоров. Далее, модельные частицы каждой из подобластей разделяются между процессорами связанной с этой подобластью группы равномерно, вне зависимости от координаты. На рис. 4 показана схема декомпозиции области для К = 4, КТ = 4.
• ■ • ■ • • ■
+
■ ф ♦ ♦ ♦
♦ ♦ ♦
+ ■ • щ +
+ 4- • ■ • • + ■ + .
♦ ♦ ■ ♦ + ■ ♦ +
• ■ +
Рис. 4. Декомпозиция области.
Область решения разбита вдоль координаты X на четыре подобласти, частицы каждой подобласти равномерно распределены между четырьмя процессорами независимо от координаты. Символы круг, квадрат, ромб, звезда означают принадлежность частиц к разным процессорам.
Ускорение вычислялось для задачи следующего размера: сетка 128 узлов по всем трем измерениям, 50 частиц в ячейке, 10 тыс. временных шагов, 12 часов счета на 32 процессорных ядрах (четырехх ядерные процессоры Xeon). Это означает, что вычисляется ускорение для 64, 128 и 256 ядер по сравнению с 32. В данном случае KL оставалось постоянным, увеличивалось только KT. Это означает, что для указанной сетки расчетная область всегда делилась на 32 подобласти, а частицы каждой подобласти разделялись между 1, 2, 4 или 8 процессорами, получено ускорение 1, 2, 3.96, 7.85 соответственно.
Вычислительные эксперименты проводились на сетках следующего размера: 128, 256 узлов по каждому измерению. Число модельных частиц в каждой ячейке от 50 до 250. Для расчетов использовались суперкомпьютеры НКС-30Т (Новосибирск, ИВМиМГ СО РАН) - до 256 ядер в отдельном расчете, СКИФ Cyberia (Томск, ТГУ) - до 500 ядер в отдельном расчете, МВС-100000 (Москва, МСЦ), до 4096 ядер в отдельном расчете.
4. РЕЗУЛЬТАТЫ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ
Для вычислительных экспериментов по моделированию взаимодействия плазмы с пучком были заданы следующие начальные условия: температура электронов 10 эВ, температура ионов 0, масса иона равна массе электрона, плотность плазмы 1017 см-3, отношение плотностей электронов пучка и электронов плазмы 10-4, энергия пучка 1 МэВ, размер области 50 RD (RD -дебаевский радиус, в данном случае RD = 7.4 • 10-6 см).
Число узлов сетки по каждому из трех изменений: MX = 128, MY = 128, MZ = 128. Число модельных частиц в каждой ячейке - 50 каждого сорта (электронов плазмы, ионов и электронов пучка). Частицы всех трех сортов равномерно распределены по расчетной области, таким образом в начальный момент задается постоянная плотность. В вычислительных экспериментах на суперЭВМ была получена модуляция плотности плазмы с амплитудой до 3 (относительно начальной плотности плазмы). Под модуляцией плотности в данном случае подразумевается возникновение в плазме с первоначально равномерной плотностью областей с повышенной либо пониженной плотностью.
ЗАКЛЮЧЕНИЕ
В работе показано, что температура при моделировании с помощью метода частиц в ячейках может быть вычислена как дисперсия импульсов модельных частиц. Исследована зависимость температуры от числа модельных частиц в ячейке, зависимость относительного изменения температуры от шага сетки. Показано, что эти зависимости являются сходящимися. Также исследована зависимость саморазогрева модельной плазмы от шага сетки. Таким образом показано, что методика вычисления температуры по модельным частицам работает неза-
висимо от шага сетки. Для вычисления температуры по предложенной методике необходимо большое количество модельных частиц и сетки большого размера. Поэтому расчеты проводились с использованием суперЭВМ.
[1] Григорьев Ю.Н., Вшивков В.А. Численные методы «частицы в ячейках». - Новосибирск: Наука, 2000.
[2] Бедсел Ч., Лэнгдон Б. Физика плазмы и математическое моделирование. - М.: Мир, 1989.
[3] Бурдаков А.В., Поступаев В.В. Особенности переноса тепла при пучковом нагреве плазмы в экспериментах на установке ГОЛ-3. Препринт № 92-9. - Новосибирск: ИЯФ СО РАН, 1992.
[4] Вшивков В.А., Вшивков К.В., Дудникова Г.И. Алгоритмы решения задачи взаимодествия лазерного импульса с плазмой // Вычислительные технологии. - 2001. - Т. 6, № 2.
[5] Кролл Н., Трайвелпис А. Основы физики плазмы. — М.: Мир, 1975.
[6] Hockney R.W. Measurements of Collision and Heating Times in a Two-Dimensional Thermal Computer Plasma // J. of Comp. Physics. - 1971. - № 8. - Р. 19-44.
Вшивков Виталий Андреевич — профессор кафедры параллельных вычислительных технологий Новосибирского государственного технического университета, д-р физ.-мат. наук
E-mail: malysh@ssd.ssec.ru Тел. 8(383) 330-96-65
Снытников Алексей Владимирович — научный сотрудник Института вычислительной математики и математической геофизики СО РАН, канд. физ.-мат. наук Тел. 8(383) 330-96-65
V.A. Vshivkov, A.V. Snytnikov
Temperature evaluation in supercomputer simulation of high-temperature plasma by the PIC method
Numerical simulation of heat conductivity in high-temperature plasma is being conducted for the 3D case. The model is based on the PIC method. Parallel implementation is built with the domain decomposition technique. The mixed Euler-Lagrangian decomposition is employed. The most difficult question is the evaluation of temperature in plasma. A technique for temperature evaluation in the PIC method is proposed, and the validity of the technique is proved.
Key words: numerical simulation, plasma instabilities, PIC method, parallel computations.
