Моделирование электромагнитных полей радиационного происхождения на многопроцессорных вычислительных системах
|
Рис 1. Равномерное по времени распределение частиц по ячейкам
Каждая частица создает
плотность тока . В узле вклад в ток дают частицы
2 и 3 с весами 2/3 (обратно пропорционально расстоянию до узла). Общий ток, таким
образом, равен . Для узла вклад дают частицы 3,
4 и 5. Частицы 3 и 5 – с весом 1/3, а частица 4 – с весом 1. Общий ток
оказывается равным . Это означает, что уравнение непрерывности выполнено интегрально,
но не выполнено ни в одной точке.
В случае равномерной
сетки проблема может быть решена следующим образом. Будем считать, что частица
представляет собой «облако» с равномерно распределенным зарядом, размер
которого равен размеру ячейки. На каждом шаге по времени это облако сдвигается.
Будем определять ток через границы ячейки, исходя их того, какой заряд прошел
через каждую границу за шаг по времени. Для представленного примера результат
будет следующий: сдвинувшись из положения 1 в положение 2, частица проносит через
границу две трети своего
заряда, так же как и частицы, двигающиеся из положения 2 в положение 3 и из
положения 3 в положение 4. Полный ток в этом случае оказывается равным . Это справедливо и для других пространственных точек.
Рассмотрим алгоритм, применимый для произвольной
декартовой сетки. Пусть – область изменения
пространственных переменных. Введем разностную сетку по переменной . Расчетная область по этой переменной разбивается на Nx
ячеек узлами , (Рис. 2). Сетка
задается следующими формулами:
;
;
Разностные сетки для переменных и вводится аналогично.
Заметим, что полуцелые точки необходимы в разностных схемах для уравнений Максвелла.
Компоненты вектора плотности электрического тока необходимо вычислять
для каждой ячейки расчетной сетки в точках, показанных на рисунке 3.
Заряд частицы распределяется между
вершинами той ячейки сетки, в которой она находится, в соответствии с
интерполяционной функцией:
,
где .
Интерполяционные функции для и строятся аналогично.
Рис 3. Размещение компонент плотности электрического тока
в ячейке разностной сетки
Пусть частица за шаг по времени перемещается из точки в точку в пределах одной
пространственной ячейки (Рис 4).
Рис 4. Начальное и конечное положения частицы
при перемещении за один временной шаг
Для каждой вершины с индексами обозначим через разность зарядов в
этой вершине в начальный и конечный моменты времени. Для непрерывности заряда потребуем,
чтобы сумма токов на ребрах, выходящих из каждой вершины (с учетом знаков),
была равна изменению заряда в этой вершине. Если
записать это условие для каждой вершины, получится 8 уравнений с 12
неизвестными (J):
; (2.1)
; (2.2)
; (2.3)
; (2.4)
; (2.5)
; (2.6)
; (2.7)
. (2.8)
Из этих уравнений независимы только семь,
так как последнее можно получить как сумму предыдущих, если учесть, что полное
изменение заряда равно 0. При вычислении токов используем дополнительные
соображения. Частица в течение всего шага по времени не меняет ни скорости, ни
направления. Естественно потребовать, чтобы ток, который она создает, был
сонаправлен скорости.
Для того чтобы доопределить систему,
необходимо потребовать, чтобы направление тока, который создает частица,
совпадало с направлением вектора с компонентами , , :
;
(2.9)
; (2.10)
; (2.11)
; (2.12)
где – скорость частицы.
Эти условия определяют недостающие 5 уравнений и позволяют вычислить токи.
Более того, вычисление токов можно организовать последовательно.
Условие (2.9) позволяет определить токи из уравнения (2.1).
Затем в уравнении (2.2) переносим в левую часть . Получаем уравнение относительно неизвестных :
.
Это уравнение решается при условии (2.10).
Аналогично решаются уравнения (2.3) при условии (2.11) и уравнение (2.4)
при условии (2.12). Из уравнений (2.5 – 2.7) выражаются оставшиеся
три неизвестных. Когда проекции скорости частицы на оси координат отрицательны,
исключение неизвестных происходит в другом порядке.
Если частица пересекает границу ячейки, то ее движение разбивается на две
части – до границы и после. На каждом из этапов движения частица находится
внутри одной пространственной ячейки и используется приведенный выше алгоритм.
Для определения плотности тока ток, создаваемый каждой частицей, делится на шаг
по времени и на площадь в
соответствии с таблицей, приведенной ниже.
Алгоритм использует для расчета токов только частицы, находящиеся внутри данной
ячейки. Поэтому не требуется вносить какие-либо коррективы в ток на границе расчетной
области или вблизи излучающих поверхностей. Алгоритм может быть использован для
реализации метода частиц на неравномерных сетках. Отсутствие разностного накопления
электрического поля, возникающего из-за разрыва токов, гарантировано тем, что
вычисления основаны на уравнении непрерывности.
Таблица – площади, соответствующие элементам тока
Элемент тока |
Площадь |
Элемент тока |
Площадь |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
3. Программная реализация для параллельных ЭВМ
Параллельная версия
программ разработана на Фортране 90 (с библиотекой передачи сообщений MPI)
и с использованием инструментальной системы параллельного программирования
Норма [10,11]. Параллельная версия может применяться на многопроцессорных
вычислительных системах с распределенной или общей памятью. Отметим, что часть
программ, разработанных для последовательной версии, включена в параллельную
версию за счет существующего интерфейса программ на языке Норма и языке
Фортран, что снизило трудоемкость разработки.
Общая схема
распараллеливания определяется информационными зависимостями между переменными,
вычисляемыми в соответствии с численными алгоритмами, и свойствами расчетной
сетки.
Численное решение
уравнений Максвелла с помощью разностных схем «крест» и ЛОС осуществляется в
декартовой системе координат на прямоугольной статической разностной сетке.
Фактически единственный
тип зависимостей в алгоритме схемы «крест» состоит в зависимости при вычислении
значения переменной в точке (i,j,k)
от значений переменных в точках (i+∆1,j+∆2,k+∆3). Наилучшие результаты при реализации алгоритмов
с такими зависимостями на статической прямоугольной сетке обычно дает модель
распараллеливания по данным. Параллельная реализация алгоритма схемы «крест»
получена при помощи системы Норма, которая позволяет получать выходные
параллельные программы на языке Фортран с использованием библиотеки передачи
сообщений MPI.
Алгоритм ЛОС имеет тот
же тип информационных зависимостей, что и алгоритм схемы «крест». Однако
стандартное распараллеливание по данным здесь не эффективно. Это связано с тем,
что ЛОС на каждом промежуточном временном шаге является неявной, а разностные
уравнения решаются методом прогонки. Предлагаемый подход к распараллеливанию
основан на конструкции алгоритма решения разностных уравнений ЛОС. Она
позволяет одновременно осуществлять прогонку для двумерного массива систем линейных
уравнений. Поэтому, помимо применения модели распараллеливания по данным,
разработан специальный эффективный способ реализации параллельной прогонки,
основанный на перераспределении данных - правых частей и матриц систем линейных
уравнений - между процессорами. Параллельная реализация метода ЛОС получена с
помощью системы Норма с включением программ на Фортране MPI,
составленных вручную.
Алгоритм вычисление
плотности электрического тока путем решения уравнений движения частиц
определяет тот же тип информационных зависимостей. Однако при этом количество
частиц в ячейках сетки может существенно различаться. Частицы неравномерно генерируются
в различных ячейках сетки. Эту неравномерность можно
оценить априорно, поскольку начальное распределение электронов предполагается
известным. В ходе вычислений частицы перемещаются между ячейками в соответствии
с уравнениями движения, распределяясь по ним также неравномерно, причем это
априорно оценить уже не возможно в силу нелинейности задачи – на движение
частиц влияет самосогласованное электромагнитное поле, которое определяется в
ходе того же расчета. Поэтому при расчете параметров частиц может
возникнуть дисбаланс загрузки процессоров. Реализовано динамическое
перераспределение между процессорами работ по расчету частиц, используя
следующие факты. Во-первых, расчет параметров данной частицы на каждом временном
слое никак не связан с расчетами параметров других частиц. Действительно, в
рассматриваемой модели электроны взаимодействуют только через самосогласованное
электромагнитное поле. На каждом временном слое, в соответствии с общей схемой
алгоритма метода частиц, при решении уравнений движения электромагнитное поле,
а, значит, и сила Лоренца рассматриваются как заданные величины. Во-вторых,
время расчета параметров частицы мало меняется при переходе с одного временного
слоя на другой.
Моделирование
самосогласованного электромагнитного поля подразумевает совместное численное
решение электродинамических уравнений Максвелла и уравнений движения частиц с
обменом данными на каждом временном слое. Параллельные программы, реализующие
отдельные этапы решения, включены в единую параллельную программу на основе следующих принципов:
распараллеливание по данным на линейку процессоров, номера которых
привязаны к узлам сетки по одной из пространственных координат;
перераспределение
данных в ЛОС;
динамическая
балансировка загрузки процессоров при расчете параметров частиц.
Всю совокупность данных
самосогласованного расчета можно разделить на динамические и статические
объекты. К первым относятся параметры крупных частиц. К статическим объектам
относятся параметры электромагнитного поля и плотности электрического тока,
которые вычисляются в ходе расчета в количестве, одинаковом для всех ячеек
сетки, и требуют одинакового времени расчета.
Статические объекты
распределены между процессорами равномерно.
Расчет динамических
объектов обычно занимает значительную часть ресурсов вычислительной системы. При этом, как отмечено выше, неравномерность
распределения частиц по ячейкам сетки делает невозможным распределение данных
«по месту», то есть так, чтобы расчет динамических объектов проводился на тех
же процессорах, где они появляются. Такое распределение приводит к простою
большинства процессоров и резко снижает эффективность распараллеливания.
Реализовано динамическое планирование вычислений, при котором на каждом
временном слое производится перераспределение частиц равномерно между
процессорами. Реализовано также пропорциональное перераспределение частиц между
процессорами. В этом случае, при перераспределении
частиц для расчета на следующем слое, учитывается время расчета каждой частицы
на предыдущем слое. Это актуально, если вычисление параметров различных частиц
требует различного числа операций.
Эффективность различных
способов планирования иллюстрирует рисунок 5. Обозначим время, затрачиваемое на один временной слой при расчете параметров частиц на процессоре с номером i, через . Заметим, что при
измерении время обменов информацией
не учитывается. Среднее время одного шага расчета на P процессорах равно Величина показывает “идеальное время” расчета. На
рисунке 5 = 30.16. Важной характеристикой является величина , показывающая
реальное время, затрачиваемое на одном шаге на расчет параметров частиц. Поскольку
все вычисления в задаче синхронизированы по временным шагам, очевидно, что общее
время расчета параметров частиц определяется величинами на отдельных шагах расчета.
Если частицы
рассчитываются в тех же процессорах, в которых они находятся (по месту), то
максимальное время счета приходится
на 2-й процессор и равно 51.5 сек. (минимальное время равно 6.48 сек. и
достигается на первом процессоре). Если частицы распределяются равномерно (по
числу), то максимальное время счета достигается на 2-м процессоре и равно 44.0 сек. (минимальное
время достигается на 9-м процессоре и равно 23.84.сек.). При пропорциональном
распределении максимальное время счета достигается
на 9-м процессоре и равно 30.36 сек. (минимальное время счета приходится на 1-й
процессор и равно 30.07 сек).
Рис 5. Сравнение
различных способов динамического
планирования вычислений
На рисунке 6
показано, как исходно неравномерное распределение динамических объектов между
процессорами, определенное на первом временном слое начальным распределением
электронов, выравнивается в процессе счета в результате пропорционального
планирования вычислений. Время расчета одного слоя становится практически одинаковым
для всех процессоров уже на третьем шаге.
Для
первого шага времена считаются равными, то есть на первом шаге используется
равномерное (по числу частиц) распределение.
На первом шаге максимальное
время счета достигается
на 2-ом процессоре и равно 44.3 сек. Выражение
= 44.3 - 30.24 = 14.06
сек показывает, насколько выполненное распределение отличается от
“идеального”. Для всех трех шагов выбрано одно среднее время, поскольку, как
отмечалось выше, эти времена от шага к шагу меняются не существенно. В данном
случае это позволяет не "перегружать" графическое представление.
На втором шаге, когда
при распределении используются времена расчета с предыдущего шага на каждом
процессоре, эта величина уменьшается: = 34.5 – 30.24 = 4.26 сек.
В этом случае максимальное
время счета достигается
на 1-ом процессоре. На третьем шаге выражение = 31.12 - 30.24 = 0.88
сек.
Отметим,
что на первом шага расчета минимальное время счета (23.65) достигается на 9-м
процессоре, на втором шаге минимальное время (28.85) достигается на 3-м процессоре,
на третьем шаге минимальное время (29.19) достигается на 6-м процессоре.
Рис 6. Выравнивание
времени счета одного шага по времени
на разных процессорах
Исходными
данными для работы алгоритма пропорционального распределения являются 1) общее количество частиц и 2) времена,
которые были затрачены на решение задачи расчета параметров
частиц на каждом из процессоров на предыдущем шаге расчета.
Помимо параллельных программ, реализующих
расчетную часть, реализован пользовательский интерфейс для задания исходных
данных и обработки результатов расчетов.
Таким образом, в состав разработанного программного
комплекса входят (Рис 7.):
-
управляющая программа, предназначенная для задания
основных параметров расчета и поддержки основных интерфейсов;
-
модуль генерации сеток, предназначенный для генерации в
расчетной области декартовых прямоугольных сеток;
-
модуль привязки объекта к сетке, предназначенный для
сопряжения геометрической модели и пространственной сетки;
-
модуль выбора выводимых результатов, предназначенный для
формирования массивов выходных данных вычислительного комплекса;
-
основная расчетная программа;
-
программа обработки результатов расчетов, предназначенная
для визуализации результатов расчета.
Рис 7. Структура программного комплекса
Ввод данных осуществляется на
персональном компьютере. Формируется набор файлов проекта, содержащих массивы
узлов разностной сетки, начальные распределения частиц и т.д. Проект по сети
передается на параллельную ЭВМ. Там осуществляется компиляция программ,
проводится расчет, формируются файлы с результатами, содержащие массивы
электромагнитных полей, плотностей токов и т.д. Результаты по сети передаются на
персональный компьютер, где осуществляется их визуализация и анализ.
4. Некоторые характеристики
проведенных расчетов
В данном разделе представлены некоторые временные характеристики,
достигнутые в нескольких расчетах методом частиц. Расчеты проводились на MBC-6000IM и MBC- 15000BM, установленных в
Межведомственном Суперкомпьютерном Центре РАН [12]. В решаемой задаче
использовалась сетка по пространству - 112 х 102 х 102. Число частиц и характер
их порождения – различный в различных
вариантах расчета. Расчет проводился на различном числе процессоров.
В
первом расчете задача решалась на MBC-6000IM, использовались 80 процессоров.
Проведено 600 шагов расчета. Максимальное число частиц на одном шаге расчета –
321691448. “Идеальное время” для 600 шагов расчета равно 2713 сек. При равномерном
распределении частиц по процессорам реальное время расчета равно 2954 сек., а
при пропорциональном распределении - 2754 сек.
Во
втором расчете задача решалась на 300 процессорах ЭВМ MBC- 15000BM. Максимальное число
частиц на одном шаге расчета – 894564000. Проведено 2000 шагов
расчета. “Идеальное время” для 2000 шагов расчета равно 9798 сек. Использовалось
пропорциональное распределение частиц по процессорам, реальное время расчета
составило 10155 сек. Ниже в таблице приведены характеристики для последнего
шага расчета (на последнем шаге расчета число частиц - максимальное).
Номер процессора |
1 |
50 |
100 |
150 |
200 |
Время расчета частиц (сек) |
11.62 |
11.65 |
11.66 |
11.56 |
11.73 |
Число частиц |
5426224 |
32175553 |
2422196 |
3134874 |
2909913 |
На последнем шаге расчета = 10.75сек., а = 11.91сек.
В третьем расчете задача
решалась на 200 процессорах ЭВМ MBC-15000BM. Максимальное число частиц на одном
шаге расчета – 997524000. Проведено 2500 шагов расчета. “Идеальное время” для 2500 шагов расчета
равно 22195 сек. Использовалось пропорциональное распределение частиц по
процессорам, реальное время расчета составило
22683 сек. Ниже в таблице приведены характеристики для последнего шага расчета.
Номер процессора |
1 |
50 |
100 |
150 |
200 |
Время расчета частиц (сек) |
18.31 |
17.55 |
16.98 |
17.65 |
18.03 |
Число частиц |
8172599 |
5195909 |
4983090 |
4665831 |
7541327 |
В
четвертом расчете задача решалась на 80 процессорах ЭВМ MBC-6000IM. Максимальное
число частиц на одном шаге расчета – 379101040. Проведено 4000 шагов расчета.
“Идеальное время” для 4000 шагов расчета равно 21794 сек. Использовалось
пропорциональное распределение частиц по процессорам, реальное время расчета
составило 22034 сек. Характеристики для последнего шага расчета:
Номер процессора |
1 |
26 |
53 |
80 |
Время расчета частиц (сек) |
10.15 |
10.17 |
10.19 |
10.19 |
Число частиц |
4439257 |
4919938 |
4967874 |
4518369 |
Для этого расчета полное
время решения задачи - 32070 сек. (включает расчет полей, расчет частиц,
накладные расходы на обмены). Полное время, затраченное на расчет частиц –
26990 сек. (с учетом времени, затраченное на организацию параллельного выполнения
расчета частиц, в том числе обмен данными между процессорами).
ЗАКЛЮЧЕНИЕ.
В данной
работе представлена трехмерная математическая модель радиационной генерации
электромагнитного поля, а также программный комплекс, разработанный для
реализации численных алгоритмов модели. Программный комплекс состоит из
параллельного вычислительного модуля, позволяющего проводить расчеты на
многопроцессорных вычислительных системах с распределеной памятью или с общей
памятью, а также пользовательского интерфейса, состоящего из препроцессора для
задания исходных данных и постпроцессора для визуализации результатов расчетов.
Параллельный
вычислительный модуль разработан на Фортране-90 с использованием библиотеки MPI и может функционировать в операционном окружении UNIX, WINDOWS.
Пользовательский интерфейс разработан на Си и предназначен для функционирования
в операционном окружении WINDOWS.
Опыт практических расчетов указывает два основных пути
развития представленной математической модели. Первый предполагает оптимизацию
параллельных программ, в первую очередь с целью более эффективного управления
памятью, расширение возможностей пользовательского интерфейса по визуализации
результатов расчетов. Второй состоит в разработке вариантов реализации для
метода частиц в конкретных задачах.
1. Д.М. Иващенко, А.А. Федоров. Российские
ускорители электронов, использующиеся в качестве моделирующих установок. – Вопросы атомной науки и техники, серия «Физика радиационного воздействия на
радиоэлектронную аппаратуру», 2002.- вып. 3, c. 120-128.
2. В.И.Бойко,
В.А.Скворцов, В.Е.Фортов, И.В.Шаманин. Взаимодействие импульсных пучков заряженных
частиц с веществом. – М.: Физматлит, 2003.
3. В. Гайтлер.
Квантовая теория излучения. – М.: ИЛ, 1956.
4. Hockney R.W., Eastwood J.W. Computer Simulation Using
Particles.– McGraw-Hill,
New York, 1981.
5. А.К. Савинский.
Взаимодействие электронов с ткане-эквивалентными средами. –
М.: Энергоатомиздат, 1984.
6. А.В.
Березин, Н.С. Келлин, М.Б. Марков, С.В. Паротькин. О столкновениях электронов в
задачах генерации электромагнитного поля потоком электронов. – Препринт
ИПМ им. М.В. Келдыша РАН, 2003, №61.
7. И. Мак-Даниель.
Процессы столкновений в ионизованных газах. – М.: МИР, 1967.
8. А.В.Березин, Н.С.Келлин, М.Б.Марков, С.В.Паротькин,
А.В.Сысенко. О математических моделях вторичной ионизации. – Препринт
ИПМ им. М.В. Келдыша РАН, 2002, №29.
9. А.В. Березин,
М.Б. Марков, Б.Д. Плющенков. Локально-одномерная разностная схема для
электродинамических задач с заданным волновым фронтом. – Препринт ИПМ
им. М.В. Келдыша РАН №31, 2005.
10. А.Н.Андрианов, А.Б.Бугеря,
К.Н.Ефимкин, И.Б.Задыхайло. Норма. Описание языка. Рабочий стандарт. Препринт
ИПМ им. М.В. Келдыша РАН, N120, 1995, 50 с.
11. www.keldysh.ru/pages/norma
12. www.jscc.ru