CC BY
62
УДК 537
Применение технологии NVIDIA CUDA
для численного моделирования распространения
электромагнитных импульсов
И.П. Молостов 2, В.В. Щербинин1
1 Алтайский государственный университет (Барнаул, Россия)
2 Институт физики им. Л.В. Киренского Сибирского отделения Российской академии наук (Красноярск, Россия)
Application of NVIDIA CUDA Technology for Numerical Simulation of Electromagnetic Pulses Propagation
I.P. Molostov1 2, V.V. Scherbinin1
1 Altai State University (Barnaul, Russia)
2 L.V. Kirensky Institute of Physics, Siberian Branch
of the Russian Academy of Sciences (Krasnoyarsk, Russia)
Рассмотрено численное решение задачи распространения электромагнитного импульса в двумерной прямоугольной области, ограниченной идеально проводящими стенками, содержащей идеально проводящий объект произвольной формы. Вычислительный алгоритм реализован на основе метода конечных разностей во временной области с использованием программно-аппаратной технологии NVIDIA CUDA. Исследованы особенности реализации метода конечных разностей для расчетов на графическом процессоре, в частности рассмотрен способ организации хранения в видеопамяти массивов значений компонент поля в точках дискретизированного пространства. Представлены численные результаты моделирования рассеяния импульса с плоским фронтом в квадратной области, содержащей цилиндрический объект. Проведен сравнительный анализ производительности разработанного алгоритма при выполнении на центральном и графическом процессорах. Показано, что для данного алгоритма применение вычислений на графических процессорах обеспечивает существенный прирост производительности.
Ключевые слова: распространение электромагнитных импульсов, вычисления общего назначения на графических процессорах.
DOI 10.14258/izvasu(2015)1.1-06
In this paper, a numerical solution of electromagnetic pulse propagation for the case of two-dimensional rectangular area bounded by perfectly conducting walls with included perfectly conducting various-shape object has been considered. A numerical algorithm has been developed using finite difference in time-domain with NVIDIA CUDA technology implementation. A finite difference method implementation for the case of computation by graphical processor has been considered. The results of numerical modeling of a plane wave pulse scattering in the square waveguide with included cylindrical object has been shown. Comparative analysis of the performance for CPU and GPU has been carried out. It is shown, that the FDTD method computation by graphical processor provides a sufficient gain in the performance.
Key words: electromagnetic pulse propagation, general-purpose graphics processing units.
1. Введение. В классической электродина-
мике расчет токов и полей в системе сводится к решению краевой математической задачи для системы дифференциальных уравнений Максвелла [1]. Полученную краевую задачу можно решать непосредственно, а можно переформулировать в ви-
де интегрального уравнения. Оба подхода имеют свои достоинства и недостатки.
Первый подход реализуется, как правило, численно. До 1980-х гг. он не имел широкого распространения, поскольку вычислительные возможности доступных компьютеров не позволяли за
приемлемое время получить решение какой-либо практически интересной задачи. Второй подход позволяет после ряда преобразований свести интегральное уравнение к алгебраическому уравнению или системе алгебраических уравнений. Это удается, если поверхности, на которых заданы краевые условия, совпадают с координатными поверхностями какой-либо ортогональной системы координат, а сами краевые условия относительно простые (например, нулевые). К сожалению, учет конечной проводимости материала или задание границ произвольной формы приводят к значительному усложнению решения и существенно затрудняют применение данного подхода во многих случаях.
Подход, использующий прямое численное решение исходной краевой задачи, в настоящее время продолжает активно развиваться. Наиболее широкое распространение в литературе получили методы конечных разностей во временной (FDTD) [2; 3] и частотной областях (FDFD) [4] и метод конечных элементов (FEM) [5]. Для решения нашей задачи нами был выбран метод FDTD.
Наибольшее распространение при решении внутренних задач электродинамики получил метод FDTD. К недостаткам этого метода можно отнести, во-первых, необходимость задания условий на границе области решения, и соответственно при решении задачи в открытой области необходимо искусственно вводить границы и задавать на них поведение полей; во-вторых, разделение области пространства на малые конечные области приводит к необходимости численного решения систем линейных алгебраических уравнений очень высокого порядка. Первый недостаток полностью устраняется при решении внутренних задач электродинамики. Что же касается второго недостатка, то он является принципиально неустранимым, однако наращивание вычислительной мощности ЭВМ и применение эффективных численных методов решения систем алгебраических уравнений позволяет обеспечить решение практически интересных задач за разумное время.
Перспективной вычислительной технологией решения систем линейных уравнений является использование вычислений общего назначения на графических процессорах (GPGPU) [6], поскольку данный класс задач хорошо поддается распараллеливанию, а графические процессоры (GPU) представляют собой системы, архитектура которых оптимизирована для параллельной обработки данных.
2. Постановка задачи. Рассмотрим бесконечный прямоугольный волновод, расположенный вдоль оси z. Вдоль оси волновода расположен цилиндр радиуса Д. Стенки волновода и цилиндра внутри него идеально проводящие. Линейный размер волновода вдоль оси x равен а, а вдоль
оси у — Ь. Волновод не имеет диэлектрического заполнения. В силу симметрии по оси г задача является двумерной. Геометрия задачи представлена на рисунке 1. Хорошо известно [1], что в такой постановке решения для компонент Ех ,Еу ,ЕХ и Нх, Ну, Нг независимы друг от друга и электромагнитное поле может быть разложено на поперечное электрическое (ТЕ) и поперечное магнитное (ТМ). Для определенности в дальнейшем будем рассматривать только ТМ поле. Возбуждение волновода происходит в начальный момент времени в окрестности оси х импульсом вида:
Ez (x,y) =
1500 если
20 <~г— < 50 и
Ay ;
x;
500 <— < 3500 Ax
(1)
иначе.
Рис. 1. Геометрия задачи
Рис. 2. Расположение различных компонент ЭМ поля в дискретизированном пространстве
0
(x, y, z, t) = (iAx, jAy, kAz, nAt),
(2)
любую функцию, зависящую от координат и времени, можно записать в виде:
F (x,y,z,t) = F n(i,j,k),
(3)
где г, к, п — целочисленные индексы. В дальнейшем будем полагать Ах = Ау = Аг. Для корректной работы метода требуется, чтобы линейные размеры элементарной ячейки были больше, чем расстояние, проходимое волной за промежуток времени Аt [7], т.е.:
УЧА^)2 + (АУ)2 + (А1)2 > Г,, ,, А/. (4)
где Стах — это максимальная скорость света в рассматриваемой области пространства. Для получения расчетных выражений следует учесть,
что в рассматриваемой задаче отсутствует зависимость поля от координаты 2, диэлектрическая е и магнитная ц проницаемости являются константами, а сторонние токи отсутствуют, = 0. Поскольку рассматривается только ТМ компонента электромагнитного поля,
Ex = Ey = Hz =0.
(5)
С учетом сделанных предположений векторные уравнения Максвелла можно расписать следующим образом:
Рис. 3. Расположение областей, в которых определены различные компоненты ЭМ поля в двумерном дискре-тизированном пространстве
3. Метод решения задачи. Метод ЕБТБ — это один из наиболее распространенных численных методов электродинамики, основанный на дискретизации уравнений Максвелла, записанных в дифференциальной форме. Одним из первых исследований, в которых использовался данный метод, была работа [7]. Согласно этому методу, непрерывно изменяющиеся в пространстве компоненты поля заменяются на значения компонент поля в узлах дискретной сетки, причем узлы для различных компонент поля смещены по отношению друг к другу на половину шага дискретизации по времени и по каждой из пространственных координат. Расположение точек, в которых необходимо вычислять компоненты поля в пространстве, представлено на рисунке 2.
3.1. Переход к конечно-разностным уравнениям. Обозначая точки дискретизи-рованного пространства-времени следующим образом:
дВх ' dt дВу dt dDz dt
dEz dy '
dEz
dx '
(6)
dHy dHx
dx
dy
Заменяя частные производные конечными разностями и применяя материальные уравнения, можно выражения (6) привести к виду:
e (i,j) = -щ-1 (i,j)+
+Z -Z
AT Ax AT ~Äy
Я" 2
Hx
i,j +
1
— Hx
1
i, j
Hxx
1 At
Ti-n + i Hy
i +
+
1 At
(7)
где At = cAt = J—At, Z = Уравнения (7)
представляют собой расчетные выражения.
3.2. Реализация разностного алгоритма
на GPU. Как видно из (7), для расчета напряженности поля в какой-либо точке дискретизиро-ванного пространства требуется только значение напряженностей полей в соседних с ней точках в предыдущий момент времени. Следовательно, для текущего момента времени поля во всех точках пространства могут быть рассчитаны независимо друг от друга, таким образом, задача легко поддается распараллеливанию.
Из рисунка 2 и выражений (7) следует, что различные компоненты электрического и магнитного полей определяются в различных точках, следовательно, можно разделить точки дискретизиро-ваного пространства на три массива: 1. Массив точек, в которых определена компонента Ez.
2
2
2
1
Рис. 4. Распределение компоненты поля Ez в пространстве при: a) t = 1000At, b) t = 1500At, c) t = 2000At, d) t = 3000At
2. Массив точек, в которых определена компонента Hx.
3. Массив точек, в которых определена компонента Hy.
При этом точек, где определена компонента Ez, меньше, чем точек, в которых определены компоненты Нх и Ну, на и соответственно. Рассмотрим пример, иллюстрирующий данное утверждение. Пусть рассматриваемая задача решается в квадратной области со сторонами 4Дх и 4Ду. На рисунке 3 показано распределение различных компонент поля в этой области пространства. Поскольку граница области является идеально проводящей, то должно выполняться условие Ez =0 для всех точек, принадлежащих ей. Следовательно, расчет компоненты Ez необходимо выполнить только в четырех внутренних точках пространства. Таким образом, точек, в которых определены компоненты Hx и Hy, по шесть на каждую компоненту.
Особенностью вычислений на GPU является необходимость при вызове функции явно указы-
вать объем данных, который будет обработан. Поскольку массивы, содержащие значения компонент поля в точках пространства, имеют различный размер, то для каждой из компонент поля необходимо использовать отдельную вычислительную процедуру.
Таблица 1
Время, требуемое для вычислений, на CPU и GPU
Количество шагов по времени, п Время выполнения на GPU, с Время выполнения на CPU, с CPU/ GPU, %
100 23.3 63.9 273.9
400 60.7 246.3 405.9
700 97.4 432.3 443.9
1500 193.3 908.3 469.3
Как уже говорилось выше, разностные методы позволяют получать решение для граничных
условий, заданных на произвольной поверхности. В случае идеально проводящего объекта, например, требуется обеспечить равенство нулю электрического поля в объеме, занимаемом этим объектом. Для реализации этого условия можно применить дополнительный массив — «маску». Каждой точке, где определена компонента Ez, соответствует одно из двух значений «маски». Если эта точка принадлежит идеально проводящему объекту, то «маска» равна 0, иначе — 1. При расчетах компоненты Ez полученный результат каждый раз умножается на «маску». В результате поле не проникает в идеальный проводник, что обеспечивает выполнение граничных условий.
4. Численные результаты. Расчетная программа была написана на специальном диалекте языка Си, предоставляемом компанией NVIDIA, реализующей поддержку CUDA API, для компиляции кода, выполняемого на графических процессорах, использовался компилятор nvcc. Код, выполняемый на центральном процессоре, подготовлен с помощью компилятора gcc. Для расчетов был использован персональный компьютер с микропроцессором Intel Core i5-2500K (3.3 ГГц) и графическим ускорителем GeForce GTS 450.
На рисунке 4 приведены графики, иллюстрирующие распределение компоненты поля Ez в пространстве в различные моменты времени. На графиках цветом изображена величина на-
пряженности компоненты поля Ez. Шкала, связывающая оттенок с напряженностью, приведена справа от каждого графика. Для наглядности рассеивающий объект выделен черным цветом, хотя напряженность на его поверхности и внутри него равна нулю.
5. Заключение. Наши результаты хорошо согласуются с результатами, полученными в работе [7]. Выигрыш в производительности, относительно расчета на одном ядре CPU, представлен в таблице 1. Видно,что реализация алгоритма FDTD с помощью CUDA API дает более высокие показатели производительности, чем реализация того же алгоритма на CPU. При распараллеливании расчетной программы на CPU, с учетом наличия у процессора i5-2500K четырех вычислительных ядер, удалось бы получить не более чем четырехкратный прирост производительности. Поскольку применяемый при расчетах GPU относится к системам начального уровня и имеет относительно низкую производительность, использование GPU представляется экономически эффективным решением для аппаратного ускорения алгоритма FDTD и других подобных численных методов.
Использование GPGPU позволяет существенно повысить производительность рабочих станций, следовательно, расширяет класс задач, решение которых не требует использования вычислительных кластеров.
Библиографический список
1. Марков Г.Т., Чаплин А.Ф. Возбуждение электромагнитных волн. — М., 1967.
2. Тихонов А.Н., Самарский А.А. Уравнения математической физики — М., 1977.
3. Поттер Д. Вычислительные методы в физике — М., 1975.
4. Albani M., Bernardi P.A. A Numerical Method Based on the Discretization of Maxwell Equations in Integral Form // Microwave Theory
and Techniques, IEEE Transactions on. — 1974. — V.22, №4.
5. Зенкевич О. Метод конечных элементов в технике — М., 1975.
6. Ong C.Y. et al. Speed it up // Microwave Magazine., IEEE. — 2010. — V.11, №2.
7. Yee K. Numerical solution of initial boundary value problem involving Maxwell's equations in isotropic media // Trans. Antennas Propagat., IEEE. — 1966. — V.14, №2.
