Доклады Российской академии наук. Математика, информатика, процессы управления, 2022, T. 505, № 1, стр. 19-23

1. ВВЕДЕНИЕ

Преимуществом гибридных суперкомпьютеров с ГПУ является высокая скорость вычислений. Она обусловлена большой пиковой производительностью ГПУ при однотипной обработке большого объема данных. Поэтому ГПУ являются перспективным средством повышения эффективности математического моделирования физических процессов [1].

Архитектура ГПУ предъявляет особые требования к вычислительным алгоритмам. Поэтому большинство научных и технических проблем, требующих численного решения уравнений математической физики, по-прежнему решаются на центральных процессорных устройствах (ЦПУ).

Данная работа представляет реализацию вычисления распределений электронов в самосогласованном электромагнитном поле методом частиц [2] на гибридном кластере. Метод подразумевает расчет электронных координат и импульсов, изменяющихся под действием поля и столкновений, совместно с решением уравнений Максвелла с плотностью тока в правой части. Гибридный кластер состоит из вычислительных узлов, в которых на каждый ГПУ приходится n штук ЦПУ. ГПУ включает m потоковых мультипроцессоров – специализированных чипов, каждый из которых является независимым вычислителем. Узел имеет общую оперативную память, доступную всем ЦПУ. Оперативная память ГПУ доступна всем его чипам. Объем вычислений динамических переменных и плотности тока превышает потребный для решения уравнений Максвелла. Поэтому использование ГПУ для реализации метода частиц является актуальной проблемой [3–8].

Данная работа представляет реализацию метода частиц на гибридном кластере для задач со следующей спецификой. Электроны генерируются во всей расчетной области и рассеиваются с малой передачей импульса. Длина пробега, ларморовский и дебаевский радиусы электронов малы по сравнению с размером расчетной области.

2. ПОСТАНОВКА ЗАДАЧИ И ЧИСЛЕННЫЙ АЛГОРИТМ

Рассмотрим кинетическое уравнение для зависящей от времени t функции распределения электронов $f = f\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ в фазовом пространстве координат r и импульсов p [9]:

где ${\mathbf{v}}$ – скорость, ${\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$ и ${\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$ – напряженности электрического и магнитного полей, ${\text{di}}{{{\text{v}}}_{{\mathbf{p}}}}$ – дивергенция в пространстве импульсов, e – заряд электрона, $с$ – скорость света, ${{\sigma }^{t}}$ и $\sigma \left( {{\mathbf{p,p}}{\kern 1pt} '} \right)$ – полное и дифференциальное макроскопические сечения рассеяния электронов, $Q\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ – интенсивность генерации электронов в фазовом пространстве; ${\mathbf{j}}\left( {t,{\mathbf{r}}} \right) = e\int {{\mathbf{v}}fd{\mathbf{p}}} $ – плотность электрического тока. Интеграл столкновений в (1) моделирует тормозное излучение, упругое рассеяние электронов, ионизацию и возбуждение среды. При низких энергиях учитывается прилипание к молекулам и рекомбинации с ионами [10–12].

Основу алгоритма решения уравнения (1) составляет представление решения суммой частиц ${{\omega }_{l}}\delta ({\mathbf{r}} - {\mathbf{r}}_{l}^{s})\delta ({\mathbf{p}} - {\mathbf{p}}_{l}^{s})$, где функции ${\mathbf{r}}\, = \,{\mathbf{r}}_{l}^{s}(t)$, ${\mathbf{p}} = {\mathbf{p}}_{l}^{s}\left( t \right)$ описывают траекторию частицы [13], появившейся в момент времени ${{t}_{l}}$ с начальными координатами $\left( {{{{\mathbf{r}}}_{l}},{{{\mathbf{p}}}_{l}}} \right)$ в фазовом пространстве. Величина ${{\omega }_{l}}$ определяет количество электронов, моделируемое частицей с номером l. Параметры ${{t}_{l}},{{{\mathbf{r}}}_{l}},{{{\mathbf{p}}}_{l}}$ и ${{\omega }_{l}}$ выбираются так, чтобы максимально точно аппроксимировать источник электронов $Q$ минимальным количеством частиц. Алгоритм сочетает статистическое моделирование рассеяния электронов с решением уравнений их движения в электромагнитном поле между столкновениями [14]. Пусть в моменты времени ${{t}_{l}} \in \left[ {{{t}_{n}};{{t}_{{n + 1}}}} \right]$ в ячейке разностной сетки образуется L частиц, в том числе, оставшихся с предыдущего временного шага с ${{t}_{l}} = {{t}_{n}}$, генерируемых источником в течение интервала $\left( {{{t}_{n}},{{t}_{{n + 1}}}} \right)$, прибывающих из соседних ячеек и образующихся при ударной ионизации. Алгоритм вычисляет функцию $f\left( {{{t}_{{n + 1}}},{\mathbf{r}},{\mathbf{p}}} \right)$, удовлетворяющую уравнению (1) с источником ${{Q}_{n}}\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ = $\sum\limits_{l = 1}^L {{{\omega }_{l}}\delta \left( {t - {{t}_{l}}} \right)\delta \left( {{\mathbf{r}} - {{{\mathbf{r}}}_{l}}} \right)\delta } \left( {{\mathbf{p}} - {{{\mathbf{p}}}_{l}}} \right)$.

Выполнение алгоритма состоит в построении для частицы траектории, которую она проходит за время ${{t}_{{n + 1}}} - {{t}_{l}}$. При этом моделируются события: выход из ячейки, поглощение, рассеяние, генерация вторичной частицы. Траектория разбивается на сегменты между событиями. Вычисляется время t_p, за которое частица проходит сегмент. Оно определяется как минимальное из интервала ${{t}_{a}} = {{t}_{{n + 1}}} - {{t}_{l}}$, оставшегося до конца шага по времени; времени до пересечения границы ячейки ${{t}_{m}}$; времени до столкновения ${{t}_{с}}$. После определения t_p интегрируются уравнения движения в заданном поле. Значения его компонент интерполируются в точку, где находится частица. Рассчитывается плотность тока, созданного частицей за время t_p. Значение ${{t}_{a}}$ пересчитывается: ${{t}_{a}} \to {{t}_{a}} - {{t}_{p}}$. Если ${{t}_{p}} = {{t}_{с}}$, то после вычисления новых координат и импульсов моделируется столкновение. Импульс частицы мгновенно изменяется статистически в соответствии с дифференциальным сечением рассеяния. Если при столкновении возникла новая частица, то она заносится в список необработанных с координатами, весом и параметром ${{t}_{a}}$ рассеянной частицы. Если после окончания построения сегмента величина ${{t}_{a}}$ остается положительной, то алгоритм расчета повторяется для следующего сегмента, начиная с определения t_p. Если ${{t}_{a}} = 0$, то частица считается достигшей верхнего временного слоя и заносится в список обработанных с текущими параметрами. Моделирование столкновений, решение уравнений движения, расчет плотности тока рассмотрены в работах [9, 14, 15].

3. РЕАЛИЗАЦИЯ АЛГОРИТМА НА ГИБРИДНОЙ ВЫЧИСЛИТЕЛЬНОЙ СИСТЕМЕ

Рассмотренный алгоритм решения уравнения (1) позволяет выделить основной объем логически простых вычислений – расчет параметров частиц, и выполнить их на ГПУ. Уравнения Максвелла (2) решаются с помощью явной разностной схемы [16] и требуют существенно меньшего объема вычислений. Поэтому задача совместного решения уравнений (1)–(2) разбивается на два этапа. Траектории частиц и плотность тока в заданном электромагнитном поле вычисляются на ГПУ под управлением ЦПУ. Затем ЦПУ вычисляют электромагнитное поле на следующем шаге. Каждое ЦПУ обрабатывает одну из прямоугольных подобластей расчетной области. Обмен данными на границах подобластей осуществляется в модели “матрица процессоров”.

Рассмотрим вычисление плотности тока. Оперативная память и чипы каждого ГПУ делятся между n ЦПУ. В части ГПУ, выделенной данному ЦПУ, обрабатываются частицы из его подобласти. На каждом временном слое ЦПУ присваивает вновь появившимся частицам параметры ${{\omega }_{l}}$, ${{{\mathbf{r}}}_{l}}$, ${{{\mathbf{p}}}_{l}}$, ${{t}_{l}}$ и передает их на ГПУ. Также на ГПУ передаются значения сеточных компонент поля на данном слое, используемые вместе с сечениями рассеяния для моделирования движения частицы всеми чипами. На одном временном шаге используемые данные не меняются. Кроме того, каждому чипу отводится свой участок общей памяти ГПУ для хранения вычисляемых данных – сеточных значений плотности тока. Это исключает синхронизацию записи вкладов частиц, обрабатываемых отдельными чипами, в вычисляемые величины. При этом объем памяти, хранящей вычисляемые данные, резко увеличивается.

Рассмотрим выполнение вычислений на ГПУ. Под вычислительным ядром понимается задание на обработку некоторой совокупности частиц. Под обработкой понимается пересчет координат и импульсов частицы при построении очередного сегмента траектории. Ядро есть совокупность исходных данных о координатах и импульсах частиц на текущем слое, применяемого алгоритма, параллельной конфигурации вычислителей. Кодом ядра называется последовательность арифметических операций над характеристиками частицы, реализующая алгоритм для каждой из них. Вычислительные потоки, выполняющие код, группируются в блоки в соответствии с конфигурацией ядра. Конфигурация задает количество вычислителей в блоке и количество блоков. Блок потоков целиком выполняется одним чипом. Чтобы потоки могли обрабатывать независимые данные, т.е. разные частицы, при выборе данных они ориентируются на номер блока и свой номер в нем. Каждый блок в комплекте вычисляемых данных хранит значения плотности тока в ячейках сетки. В процессе обработки потоки блока атомарно записывают вклады частиц в вычисляемые данные для соответствующих ячеек. Вклады частиц блоков по окончании моделирования на данном временном слое суммируются.

За один запуск ядра для частицы строится очередной сегмент траектории. Для отдельной частицы результатом запуска может быть: переход на следующий слой и окончание моделирования на данном шаге по времени; достижение границы ячейки или точки взаимодействия; поглощение и окончание моделирования. Количество обрабатываемых при однократном запуске ядра частиц есть произведение количеств потоков в блоке N_block = 512 и числа мультипроцессоров N_mp = m/n, т.е. несколько тысяч. Для перевода всех частиц на следующий слой ядро запускается многократно. Вычислители, частицы которых после очередного запуска обработаны, при следующем запуске начинают обрабатывать другие частицы системы. Запуски ядра продолжаются до достижения всеми частицами данного ЦПУ следующего слоя.

Данные о частицах на ГПУ хранятся в трех доменах памяти. В первом домене хранятся данные частиц, оставшихся в подобласти с предыдущего шага по времени. Для каждого из потоков закладывается собственный домен под частицы, всего N_block*N_mp доменов. Второй и третий домены хранят данные входящих и выходящих частиц. К входящим относятся частицы, рожденные на текущем шаге по времени, а также прилетевшие из соседних подобластей сетки. К выходящим относятся частицы, вылетевшие за границу подобласти. Для входящих и выходящих частиц создается по отдельному домену – “корзине” на каждый чип, всего N_mp “корзин” каждого типа. Входные корзины равномерно заполняются на стороне ЦПУ в перерывах между запусками ядра. Аналогично на стороне ЦПУ опустошаются выходные корзины. Для каждого потока на ГПУ хранятся номер I_write элемента памяти, из которого берутся характеристики частицы для обработки при запуске ядра, и номер I_read элемента памяти, в который будут записаны изменившиеся характеристики обработанной частицы. Для каждого чипа на ГПУ также хранятся количества частиц во выходной (I_out) и входной (I_in) корзинах.

Последовательность вычислений (рис. 1) исключает фрагментацию памяти в массиве, т.е. возникновения пустых элементов. Если частица достигает конца шага по времени, оставаясь в текущей подобласти, то указатели I_read и I_write одновременно перемещаются на следующий элемент памяти. Если частица выбывает из рассмотрения и в массиве и образуется свободная ячейка, то вперед двигается только указатель I_read. Таким образом, все последующие частицы после обработки сдвинутся в массиве данных и закроют “дыру”. Это позволяет автоматически балансировать количество частиц между потоками. Поток забирает на обработку частицы из входных корзин только после того, как все частицы, оказавшиеся в его домене на начало шага по времени, будут обработаны. Поэтому чем больше у потока частиц для обработки, тем позднее он начнет брать дополнительные частицы из входной корзины. Наполнение корзин балансируется на ЦПУ. После окончания обработки всех частиц на всех ГПУ вклады в плотности токов от всех блоков суммируются на каждом из ГПУ. Вклады от каждой подобласти отправляются на ЦПУ для решения уравнений Максвелла на верхнем слое. Проблемой данной реализации является необходимость хранения в памяти ГПУ отдельного комплекта вычисляемых данных для каждого чипа. Уменьшить объем памяти, занимаемой комплектом, позволяет управление каждым ГПУ не одним, а несколькими (n) ЦПУ. Пусть ГПУ обрабатывает частицы в N ячейках сетки. Каждый ЦПУ управляет обработкой частиц на ячейках сетки в количестве N/n. На один ЦПУ приходится 1/n памяти и m/n чипов ГПУ. Тогда расход памяти ГПУ на хранение вычисляемых данных пропорционален (Nm/n²).

Следующей проблемой является межпроцессорная синхронизация окончания шага по времени. Она необходима, поскольку частицы могут перемещаться между подобластями в рамках одного шага. При моделировании очередного шага на каждом ЦПУ одновременно происходят два параллельных процесса: запуск на ГПУ ядер с обменом частицами и координация MPI-обменов. При координации MPI-обменов итеративно проверяется готовность ЦПУ. Считается, что оно закончило моделирование, если все частицы его подобласти достигли конца шага по времени, а для всех частиц, отправленных в соседние подобласти, получено подтверждение о получении адресатом. В случае готовности ЦПУ выставляет асинхронный барьер, не блокирующий координацию MPI-обменов. Таким образом ЦПУ, закончивший моделирование, при получении новых частиц из соседней подобласти возобновляет их обработку. По достижении всеми ЦПУ барьера моделирование приостанавливается. Подсчитывается общее количество необработанных частиц. Если оно равно нулю, шаг по времени заканчивается.

Эффективность реализации решения уравнения (1) проверена в ИПМ им. М.В. Келдыша РАН на одном узле вычислительного кластера К60. Он включает 32 ядра ЦПУ в составе двух процессоров Intel Xeon Gold 6142 v4 и четыре ГПУ nVidia V100 GV100GL, 768 Gb RAM, 2 Тб HDD. Рассмотрена задача моделирования движения в постоянном магнитном поле электронов, образуемых однородным объемным источником. Параметры подобраны так, что частицам приходится многократно пересекать границы подобластей. Сравнивались скорости расчета параметров 30 млн. частиц только на ЦПУ и в рассмотренной модели распараллеливания с использованием ГПУ. Показано, что применение ГПУ обеспечивает ускорение более, чем в 30 раз. С уменьшением числа частиц ускорение растет. Это связано со снижением затрат ресурсов на пересылку параметров частиц между подобластями. Также уменьшение числа частиц сохраняет равномерность их распределения в пространстве, т.е. по вычислителям ГПУ.

Авторы надеются, что предложенное решение проблемы синхронизации записи результатов и обмена данными при расчете параметров частиц поможет заинтересованным специалистам при решении практических задач моделирования объемной генерации, распространения и рассеяния релятивистских электронов в самосогласованном электромагнитном поле.