Реализация разностного решения уравнений Максвелла на графическом процессоре методом пирамид Текст научной статьи по специальности «Математика»

РЕАЛИЗАЦИЯ РАЗНОСТНОГО РЕШЕНИЯ УРАВНЕНИИ МАКСВЕЛЛА НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ МЕТОДОМ ПИРАМИД

С.А. Малышева 12, Д.Л. Головашкин 1,2 1 Институт систем обработки изображений РАН, Самара, Россия, 2 Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет) (СГАУ), Самара, Россия

Аннотация

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

Ключевые слова: FDTD, уравнения Максвелла, метод пирамид, GPU, CUDA.

Цитирование: Малышева, С.А. Реализация разностного решения уравнений Максвелла на графических процессорах методом пирамид / С.А. Малышева, Д.Л. Головашкин // Компьютерная оптика. - 2016. - Т. 40, № 2. - С. 179-187. - DOI: 10.18287/2412-6179-2016-40-2179-187.

Введение

Метод разностного решения уравнений Максвелла во временной области (FDTD) [1] широко используется в современных исследованиях при моделировании: распространения излучения [2], материалов с новыми свойствами [3], взаимодействия излучения с биологическими объектами [4] и во многих других приложениях.

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

Сталкиваясь с известными трудностями («кремниевый тупик»), развитие аппаратной базы вычислительной техники в последнее десятилетие связывается не с наращиванием тактовой частоты процессоров, а со специализированными архитектурными решениями, позволяющими значительно увеличить быстродействие операций определенного типа. Наиболее заметные успехи в этом направлении достигнуты при развитии архитектуры графических процессоров (GPU) и соответствующих программных инструментов (CUDA [5], OpenCL [6]), позволяющих многократно ускорить действия с векторами и матрицами.

Современные реализации FDTD-метода на GPU [7, 8], в основе которых лежит операция вычитания матриц, характеризуются повышением производительности на порядок по сравнению с вычислениями на центральном процессоре. Однако ограниченный объем видеопамяти (2-6 Гб. по сравнению с 4-32 Гб. оперативной памяти) не позволяет в полной мере использовать преимущество GPU в быстродействии.

Снижение влияния упомянутого фактора представляется авторам весьма актуальной задачей. Для ее решения в настоящей работе предлагается развить метод пирамид [9, 10, 11], основанный на уменьшении интенсивности коммуникаций между оперативной и видеопамятью за счёт дублирования арифметических операций (в частности, ГБТБ-метода).

1. Разностная схема Уев

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

Табл. 1. Коэффициенты сеточных проекций электрического и магнитного полей

A d Mt a Ix b Jy c Kz

Ex 0,5 1

Ev 0,5 1

Ez 0,5 1

Hx 0,5 1 0,5 1 0,5 1

Hy 0,5 1 0,5 1 0,5 1

Hz 0,5 1 0,5 1 0,5 1

Рассмотрим трехмерную по пространству область вычислительного эксперимента Б3 (0 < г < Т, 0 < х < Ьх,

0 <у <Ьу, 0 < г < Ьг) с наложенной на нее сеточной структурой Б1, в узлах которой {(Ша, х++а, У]+ъ, хк+с): Ш+й = (ш+й)к, ш = 0, 1,...,. М-М, (М= Т/кг), Хг+а = (1+а)кх,

1 = 0, 1,., 1-1х, (I=Ьх/кх), у+ъ = а+Ъ)ку, ] = 0, 1,...,—у, (3 =Ьу/ку), 2к+с = (к+с)кг, к = 0, 1,...,К-КЪ К=Ь/кг] определены сеточные проекции поля АШ^й+ъ к+с.

Значения коэффициентов а, Ъ, с, й для различных сеточных проекций электрического и магнитного полей

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

Разностная схема для этой области представлена уравнениями (1)-(6):

и:

= А„

и:-

),/+0.5 ,¿+0.5 ^¡,у+0.5,к+0.5

-А

Ет — Ет Ет — Ет

г1, _/+1,к+0.5 г1, /,¿+0.5 у", /+0.5,к+1 У1, j+0.5,¿

А„

Ит+0.5 = а ит—05 — А

у+0.5,/",¿+0.5 "¡+0.5,/ ,¿+0.5 у+0.5,/ ,¿+0.5 Ь

У

-Е

г ет х

л1+0.5,/,¿+1 х1+0.5, /,к

Е"" — Ет

z¡+1, /,¿+0.5 г, / ,¿+0.5

А.

ит

=А

И:—05 — А

2Г + 1

Ет — Ет Ет — Ет

у+1,/+0.5,¿ у",/+0.5^ Х1+0.5,/+1^ Х1+0.5, / ,¿

"¡+0.5,/+0.5^ г+0.5,/+0^

"¡+0.5, / Л х+0.5,/^

и т+0.5 — ит+0.5

^¡+0.5,/+0.5,¿ гй-0.5,/—0.5^ у+0.5,/ ,¿+0.5 у+0.5,/,¿—0.5

ит

А

—и

У / т+0.5

А

Ет

у,/+0.5 ^

= с е:

+с

Ь,/+0.5 ,¿ у,/+0.5 ,¿ Ьi,/+0.5^

",/,¿+0.5 21,/ ,¿+0.5

и т+0.5 — и т+0.5

х1,/+0.5 ,¿+0.5 х1,/+0.5,¿—0.5

и т+0.5 — и т+0.5

у+0.5, /,¿+0.5 y¡—0.5, /,¿+0.5

ит

—и

т+0.5

ит

—ит

К

(1) (2)

(3)

(4)

(5)

(6)

где

Са = (1 — А/А)/а

2Ч/Л £0

К

2Ч/* е0

а ]А

¡^ £,-,м £<> / 2е,-,м е0

а * А / а * А

А = (1 —

а(,„ = (>л.+/),

/ 2ти,кт0

Табл. 2. Диапазоны индексов для проекций электрического и магнитного полей

£0, тс - электрическая и магнитная постоянные, £/

тк - сеточные относительные электрическая и магнитная проницаемости изотропной линейной бездисперсной среды, а/ а*/ - сеточные удельные электрическая и магнитная проводимости.

При моделировании бесконечного свободного пространства вокруг А3 наложим на сеточную область поглощающие слои СРМЬ [1] шириной пхРМЬ_\ отсчётов слева и пхРМЬ_2 справа вдоль оси х, пуРМЬ_1 и пуРМЬ _2, пгРМЬ_\ и пгРМЬ_2 по осям у и г соответственно. В табл. 2 представлены диапазоны индексов (в нотации Голуба [13]) для каждой сеточной функции, соответствующие подобластям с наложенными слоями. Выражения (7)-(12) описывают сеточные уравнения в поглощающих слоях.

А

¡,/+0.5^+0.5

¡,/+0.5^+0.5

г

А

х

А

¡+0.5, /+0^

х

Ет+1 = С Ет + С

ь

х

А

(+0.5, / ,¿

¡+0.5,/ ,¿

у

(+0.5,/+0.5^

(—0.5,/+0.5^

А

А

г

х

/+0.5^+0.5

(,/—0.5,¿+0.5

Ет+1 = С Е + С

ь

г

¡, / ,¿+0.5

¡, / ,¿+0.5

А

х

\ 1 ] k

Ех - - 1: пуРМЬ 1 ^пуРМЬ 2+1:7 1 :пгРМЬ 1 К-пгРМЬ 2:К

Еу 1 :пхРМЬ 1 1-пхРМЬ 2+1:/ - - 1 :пгРМЬ 1 К-пгРМЬ 2+1:К

Ег 1 :пхРМЬ 1 1-пхРМЬ 2+1:/ 1: пуРМЬ 1 7-пуРМЬ 2+1:7 - -

Их - - 1 :пуРМЬ 1-1 7-пуРМЬ 2+1:7-1 1:пгРМЬ 1-1 К-пгРМЬ 2+1:К-1

Иу 1:пхРМЬ 1-1 /-пхРМЬ 2-1:1-1 - - 1 :пгРМЬ 1-1 К-пгРМЬ 2+1:К-1

Иг 1:пхРМЬ 1-1 /-пхРМЬ 2+1 :/-1 1:пуРМЬ_1-1 7-пуРМЬ_2+1:7-1 - -

и:+05 = Ба и:—05 — аь

x¡,/+0.5,¿+0.5 ",/+0.5 ,¿+0.5 x¡,/+0.5 ,¿+0.5 ь¡,/+0.5^+0.5

и",+05 = А ит—05 — А.

y¡+0.5, /,¿+0.5 "¡+0.5, /,¿+0.5 у+0.5, /,¿+0.5 ь+0.5, /,¿+0.5

и",+05 = А ит—05 — А,

^+0.5,/+0.5 ,¿ "¡+0.5,/+0.5 ,¿ ^¡+0.5,/"+0^ ь+0.5,/"+0^

( е: — Ет Ет — Ет

Z¡,/+1 ,¿+0.5 z¡, /,¿+0.5 y¡,/+0.5 ,¿+1 y¡,/+0.5 ,¿

КУ/+05К

к А.

+ — ^:

ИХУ1, /+0.5^+0.5 Иxz¡,/+0.5,¿+0.5

Ет — Ет Ет — Ет

•х+0.5,/,¿+1 х1+0.5,/,¿

КА

-1, / ,¿+0.5 ^,/,¿+0.5 + ^^: ^^:

ИУz¡+0.5,/ ,¿+0.5 ИXУx¡+0.5,/,¿+0.5

КА

Ет - Ет Ет - Ет

у+1,/+0.5 ,¿ У¡,/+0.5^ x¡+0.5,/+1,¿ x¡+0.5,/,¿

КА

КУ/+05 АУ

,/,¿ + — ^:

0.5,/+0.5^ ИZУ¡+0.5,/+0.5 ,¿

е:+1 = с Е: + С

т

"¡+0 5 / k X( +

(ит+0.5 — ит+0 5 ит+0 5 — и:+0'5

^+0.5,/+0.5,¿ ^¡+0.5,/—0.5,¿ y¡+0.5, /,¿+0.5 y¡+ 0.5, / ,¿—0.5

'¡+0.5,/ ,¿ ¡+0.5 ,/,¿ ¡+0.5,П+0.5,/,¿

К У/АУ

КА

т+0.5 — ^ т+0.5

Ex,У¡+0.5,/,¿ Ex,z¡+0.5,/,¿

(7)

(8)

(9)

(10)

¿+0.5

¿+0.5

¡+0.5

¡+0.5

Em = Ca Em + Cb

y,, /'+0.5ai,j+ 0.5,k yi,j+0.5,k bi,_,

^ Hm+0 5 — Hm+0 5

xi, j+0.5,k+0.5 xi, j+0.5,k—0.5

к A

^ Hm+0 5 _Hm+0 5

yi+0.5,j,k 0.5 yi—0.5,j ,k+0.5

rm+0.5

Hm+05 — H

zi+0.5,j+0.5,k zi—0.5, j+0.5,k , \T/m+0.5 \T/m+0.5

k К

+ y m+05 — Y m

Vi, j+0.5,k Vi,j+0.5 ,k

i, j,k+0.5 zi, j ,k+0.5 bi, j,k+0.5

m+0.5 x,

k К

Hm+0.5 — H

xi, j+0.5,k+0.5_xi, j—0.5,k+0.5 + Y m+0.5 — Y m+0.5

к h Bzxi,j,k+0.5 Ezy,,j,k+0.5 ^

(11)

(12)

где

Ym

= b Y

yj+°5 xyi, j+0.5,k+0.5

+ c,,

i Em

zi, j+1,k+0.5

— E

Y

m—0.5

xy,+0.5,j,i JJ xyi+0.5, j,k.

+c

У / m+0.5

Hm+0.5 — H

zi+0.5, j+0.5 ,k zi+0.5, j—0.5,k

_0»_+f_w |ht

b = g l_>*». _01 w

w . w

'0 lsw c0

— 1

sw, Kw, aw - коэффициенты CFS (Complex-Frequency-Shifted)-тензора. YH , Ym+05

для остальных ком-

понент определяются аналогично с точностью до замены х, у, г и I,], к. Коэффициенты и с„ отличны от 0 только в поглощающих слоях.

По краям сеточной области зададим граничное условие Дирихле, соответствующее идеальному проводнику, когда приведенные в табл. 3 компоненты вектора электрического поля АШ+ау+Ък+с равны 0 на соответствующих границах при 0 < т < М:

Табл. 3. Граничное условие Дирихле для проекций электрического и магнитного полей

.................................. 0 < i < I-Ix 0 < j < J-Jy 0 < k < K-Kz

0 < i < I-Ix - E"; = 0, Em = 0 xi+0.5, j,0 xi+0.5, j ,K e; = 0, E; = 0

0 < j < J-Jy Em = 0 , Em = 0 yi,j+0.5,0 Ун,j+0.5,K - Em = 0, Em = 0 y0, j+0.5,k yI, j+0.5,k

0 < k < K-Kz Em = 0, Em = 0 zi ,0,k+0.5 ' z,,J ,k+0.5 Em = 0 , Em = 0 z0,j,k+0.5 У ZI, j,k+0.5 -

Em 1 = C Em + C

z

i, j ,k+0.5

m—1

i, j,k+0.5

+0.5

Vi, j+0.5,k+0.5

h

y

s

w

cw =

e

2

Sw Kw + Kwaw

Начальное условие зададим как Е0 = Е° = Е0 = 0 , что соответствует от-

х1+0.5, у,к у1,;+0.5,к Г1,;,к+0.5

сутствию излучения в начальный момент времени. Поле вводится в область вычислительного эксперимента с помощью метода результирующего поля [14].

2. Метод пирамид

Приступая к синтезу алгоритма для организации вычислений по (1-12) на графическом процессоре, остановимся на методе пирамид, позволяющем снять остроту известного ограничения на объем видеопамяти.

Суть классического метода пирамид [10] заключается в следующих предписаниях:

• построении пространства итераций циклической конструкции, обеспечивающей вычисления по (112) на всей сеточной области Е)\ ;

• выделении на этом пространстве результирующих итераций, от которых не зависит никакая другая;

• отнесении к ведению каждой задачи искомого алгоритма (в терминах модели канал/задача) всех итераций, от которых зависит данная результирующая. Проблема ограничения памяти графического вычислительного устройства решается в [11] последовательной реализацией таких задач на одной видеокарте, а именно:

• пространство итераций разбивается по переменной, отвечающей за шаги по времени, на равные интервалы, к каждому из которых применяется метод пирамид;

• необходимо, чтобы все итерации, соответствующие меньшему значению переменной, не зависели от итераций, соответствующих большему значению переменной;

• высота пирамид выбирается из соображений минимизации длительности вычислений;

• разбиение на задачи может быть одинаковым для всех проходов;

• задачи выполняются квазипараллельно, по очереди занимая вычислительное устройство на время, необходимое для выполнения одного прохода;

• каждая из задач решается в векторизованном виде. Применение модифицированного подхода к сеточным уравнениям позволяет сократить число коммуникаций, т.к. в нем пересылка данных осуществляется только после расчёта нескольких временных слоёв (их число равно высоте пирамиды), а не на каждом шаге, как в классическом алгоритме. Квазипараллельное исполнение задач позволяет одной задаче монопольно занимать GPU, наиболее полно используя доступную видеопамять.

Особенностью FDTD-метода является выражение каждой сеточной функции через значение сеточных функций на предыдущих временных слоях в явном виде, поэтому метод пирамид применяется ко всем уравнениям (1)-(6) одновременно: на каждом временном шаге рассчитываются все 6 компонент А'+^+ь k+c

электромагнитного поля; перед началом расчёта и после завершения обработки пирамиды осуществляется пересылка необходимых значений для всех компонент полей. При этом расчёт полей в поглощающих слоях по формулам (7)-(12) будет производиться только для пирамид, расположенных по краям исходной сеточной области. Количество таких пирамид определяется соотношением ширины поглощающего слоя, ширины основания пирамиды и высоты пирамиды, определяющей ширину области перекрытия.

2.1. Случай одномерной сеточной области

Проиллюстрируем применение метода пирамид к реализации схемы Yee на GPU на примере одномерной задачи. Для одномерного случая (плоская электромагнитная волна распространяется по направлению Z) от нуля будут отличны только компоненты Ex, Hy электромагнитного поля и схема Yee сократится до двух уравнений (2), (4), в которых сеточные функции других компонент обнуляются. Расчёт по этим уравнениям требует наличия в оперативной памяти всей сеточной области.

Для обхода этого ограничения в [11] предлагается разбивать сеточную область на подобласти, каждая из которых может быть помещена в память GPU и производить в ней итерации по времени стандартным образом, после чего перемещать обратно на CPU. В случае тривиального алгоритма производится одна итерация по времени, при этом для расчёта r значений Em и Hm+°;5 необходимо

передать по R = r+2 значений каждой из компонент, так как дифференциальный шаблон затрагивает соседние узлы. Однако этот алгоритм предполагает пересылку данных на каждом временном шаге.

Чтобы этого избежать, в [11] предлагается следующая модификация: для вычисления M слоёв сеточных уравнений (2), (4) необходимо последовательно выполнить s проходов методом пирамид, вычисляя каждый раз n слоёв (M = sn, n - высота пирамиды). Взаимодействие между задачами происходит через CPU и осуществляется только по окончании прохода, т.е. проходы выполняются независимо каждой из задач. При этом для расчёта в результате прохода r

значений Ex и H y необходимо передать по

хк ук

R = r + 2n значений каждой из компонент. Число задач m выбирается исходя из соотношения 2rm = 2K -1, где K - количество узлов сеточной области по направлению Z.

Множество итераций, выполняемых каждой задачей, представляет собой усеченную пирамиду высотой n, шириной верхнего основания r и нижнего основания R. Это иллюстрирует рис. 1 для второй задачи на примере разбиения области на 3 подзадачи, n = 3, r = 3, R = 9, где квадраты обозначают значения сеточной функции, столбцы - пространственные узлы сетки, строки - временные слои, цифры - номер задачи, вычисляющей значение в данном узле. На каждом временном слое вычисляются значения E и H .

xk yk+0.5

Темно-серым обозначены значения, вычисляемые второй задачей, бледно-серым - передаваемые для ее расчёта.

Оценка сверху числа значений сеточной функции, обрабатываемых за один проход одной задачей, для метода пирамид и тривиального алгоритма приведена в табл. 4, при этом для расчёта n временных слоёв тривиальным алгоритмом необходимо выполнить n проходов.

задача 1 задача 2 задача 3

npoxodl п 2

проходу п

Рис. 1. Декомпозиция сеточной области в одномерном случае для второй задачи Табл. 4. Оценка числа значений сеточной функции для одной задачи за один проход в одномерном случае

............................................................................. Тривиальный алгоритм Метод пирамид

Пересылаемых всего: 4(R-1) 4(R-n)

- в видеопамять 2R 2R

- из видеопамяти 2(R-2) 2(R-2n)

Вычисляемых по шаблону 2(R-2) 2n(R-n-1)

Таким образом, для одной задачи число пересылок сокращается в n (R -1)/(R - n) раз, а вычислительная сложность увеличивается в (R - n - 1)/(R - 2) раза.

2.2. Случай двумерной сеточной области

Иллюстрируя применение метода пирамид к двумерной сеточной области, рассмотрим пример распространения TM-волны. В этом случае от нуля будут отличны только компоненты Ex, Hy, Hz электромагнитного поля и схема Yee сократится до трёх уравнений (2), (3), (4), в которых сеточные функции других компонент обнуляются.

Тривиальный алгоритм для двумерной по пространству задачи с одномерной декомпозицией сеточной области аналогичен изложенному в п. 2.1. Сеточная область разбивается на равные прямоугольные подобласти размером J х (R - 2), которые могут быть целиком помещены в видеопамять; т.к. для вычисления одного значения требуются значения сеточной функции в соседних узлах. Подобласти должны перекрываться на два сеточных узла по каждой стороне, не прилегающей к краям сеточной области. На каждом шаге по времени производится последовательная обработка подобластей: значения сеточных функций копируются в память GPU, вычисляются значения для следующего временного слоя, вычисленные значения копируются на CPU.

Аналогично для организации вычислений по методу пирамид с одномерной декомпозицией требуется область видеопамяти для размещения по R = r+2n значений каждой из компонент электромагнитного поля Ex, Hy, Hz для вычисления по r значений для каждой компоненты. Число задач m выбирается исходя из соотношения 3rm = 3K - 1.

На рис. 2 показана одномерная декомпозиция двумерной области на примере п = 3, г = 3, К = 9, где кубы обозначают значения сеточной функции, столбцы -пространственные узлы сетки, строки - временные слои, цифры - номер задачи, вычисляющей значение в данном узле. На каждом временном слое вычисляются значения Ех.р Ну.^ Я^к .

задача 1 задача 2 задача 3

8

-Т5

¡С

1,2 1,2 2 2,3 2,3 1,2 1,2 1,2 1,2,3 2,3 2,3 2,3

1,2 1,2 2 2,3 2,3 1,2 1,2 1,2 1,2,3 2,3 2,3 2,3

R

Рис. 2. Декомпозиция сеточной области в одномерном случае для второй задачи

Темно-серым обозначены значения, вычисляемые второй задачей, бледно-серым - получаемые ею от других задач.

Оценка сверху числа значений сеточной функции, обрабатываемых за один проход одной задачей, для метода пирамид и тривиального алгоритма приведена в табл. 5, при этом для расчёта п временных слоёв тривиальным алгоритмом необходимо выполнить п проходов.

Как видно из табл. 5, количество операций отличается от одномерного случая только множителем (3/- 1), соответствующим количеству узлов по направлению у для компонент Ех, Ну, Я Таким образом, для одной задачи число пересылок, так же как и в случае одномерной задачи, сокращается в п (К -1)/(К - п) раз, а вычислительная сложность увеличивается в (К - п -1)/(К - 2) раза.

Табл. 5. Оценка числа значений сеточной функции

для одной задачи за один проход в двумерном случае

............................................................................. Тривиальный алгоритм Метод пирамид

Пересылаемых всего: 2(3J-1)(R-1) 2(3J-1)(R-n)

- в видеопамять (3J-1)R (3J-1)R

- из видеопамяти (3J-1)(R-2) (3J-1)(R-2n)

Вычисляемых по шаблону (3J-1)(R-2) (3J-1)n(R-n-1)

2.3. Случай трёхмерной сеточной области с одномерной декомпозицией

В случае трёхмерной сеточной области от нуля отличны все компоненты электромагнитного поля и схема Уее представлена уравнениями (1)-(6).

Методики построения тривиального алгоритма и алгоритма метода пирамид, применяемые к трёхмерной сеточной области при одномерной декомпозиции последней, аналогичны рассмотренным выше с той лишь разницей, что сеточные подобласти, которые выделяются задачам, имеют форму параллелепипедов размером I х /х (К - 2), которые могут быть целиком помещены в видеопамять.

Оценка сверху числа значений сеточной функции, обрабатываемых за один проход одной задачей, для метода пирамид и тривиального алгоритма приведена в табл. 3, при этом для расчёта п временных слоёв тривиальным алгоритмом необходимо выполнить п проходов.

Как видно из табл. 6, количество операций отличается от одномерного случая только множителем (61/31 - 3/ + 1), соответствующим количеству узлов в плоскости ху для компонент электромагнитного поля.

Табл. 6. Оценка числа значений сеточной функции для одной задачи за один проход в трехмерном случае

.............................................................................. Тривиальный алгоритм Метод пирамид

Пересылаемых всего юс 2(6IJ-3I-3J+1) (R-1) 2(6IJ-3I-3J+1) (R-n)

- в видеопамять (6IJ-3I-3J+1) R (6IJ-3I-3J+1)R

- из видеопамяти (6IJ-3I-3J+1) (R-2) (6IJ-3I-3J+1) (R-2n)

Вычисляемых по шаблону Юп (6IJ-3I-3J+1) (R-2) (6IJ-3I-3J+1) n(R-n-1)

На рис. 3 представлена блок-схема алгоритма расчёта по методу пирамид, вычисления по которой описаны в следующем параграфе.

3. Экспериментальное исследование метода пирамид

Ниже описаны эксперименты, проведенные с целью выявления ускорения вычислений при применении метода пирамид к FDTD.

Вычислительные эксперименты по определению ускорения проводились на видеокарте NVIDIA GeForce GTX 660 Ti (табл. 7) и процессоре Intel Core 2 Duo E8500 (табл. 8). Исследование велось в операционной системе Ubuntu.

3.1. Программная реализация FDTD-метода

Рассмотрим распространение плоской волны в трехмерной области свободного пространства, ограниченной поглощающими слоями CPML, для чего будем рассматривать уравнения (1)-(6), решение которых соответствует моделированию распространения поля в подобласти без PML; (7)-(12) - в подобласти с PML и (13) - как начальное условие.

Зададим сеточную область размером 256*256*256 отсчётов по пространству, ширину поглощающего слоя выберем nxPML = nyPML = nzPML = 11 отсчётов. Для размещения такой задачи в памяти видеокарты необходимо 467 Мб.

J

Рис. 3. Схема алгоритма метода пирамид Табл. 7. Основные характеристики GPU NVIDIA GeForce GTX 660 Ti

Характеристика Значение

Количество мультипроцессоров, шт. 7

Размер видеопамяти, Гб 2

Максимальное число потоков в блоке, шт. 1024

Максимальная размерность блока потоков (х, у, г), шт. 1024x1024x128

Максимальная размерность сетки блоков, шт. 2147483x65535x65535

Тактовая частота ядра, МГц 1032

Тактовая частота памяти, МГц 6008

Табл. 8. Основные характеристики СРиIntel Core 2 Duo E8500

Характеристика Значение

Тактовая частота ядра, ГГц 3,16

Тактовая частота шины СРи, МГц 1333

Кеш Ь1, Кб 64x2

Кеш Ь2, Кб 6144

Применим метод пирамид с одномерной декомпозицией и проведем разбиение исходной сеточной области по направлению К. Ниже приведен исходный код, описывающий алгоритм, приведённый на рис. 3:

Г

host_ void raschetGPU_pyramid_1d()

//Число запускаемых блоков

(Imax-(Jmax-

/

int SIZEx = ((Imax-1)%BLOCK_DIMx==0) 1)/BLOCK_DIMx : (Imax-1)/BLOCK_DIMx+1; int SIZEy = ((Jmax-1)%BLOCK_DIMy==0) 1)/BLOCK_DIMy : (Jmax-1)/BLOCK_DIMy+1; dim3 gridSize = dim3(SIZEx,SIZEy, 1); dim3 blockSize = dim3(BLOCK_DIMx,BLOCK_DIMy, BLOCK_DIMz);

create_temp_fields();//Копирование пересекаю-//щихся областей во временный буфер int countPyramidsK = ceil((Kmax-1) (pyramidBaseLengthK + 0.0f)); // Число пирамид (п. 2 на рисунке 3) int currentTime = 1;//Текущий временной шаг int timeOfPass =((nmax)%PASS==0)? nmax/PASS : nmax/PASS+1;

// Число проходов (п. 3 на рисунке 3) //Цикл по проходам (п. 4 на рисунке 3) for (int pass = 1; pass <= PASS; pass++)

int currentBaseLengthK; //ширина по К верх-

//него основания текущей пирамиды

int leftOffsetK; //ширина левого перекрытия

//по К нижнего основания пирамиды

int rightOffsetK; //ширина правого перекры-

//тия по К нижнего основания пирамиды

int fullBaseLengthK; //ширина по К нижнего

//основания текущей пирамиды

int leftPyramidBorderK; //индекс левой грани-//цы основания пирамиды

int rightPyramidBorderK;//индекс правой гра-//ницы основания пирамиды

int durationOfTimePass = timeOfPass * pass;

//временной шаг, соответствующий окончанию //текущего прохода

//Копируем во временный буфер начальные зна-//чения полей для прохода copyFields(...);

//Цикл по пирамидам (п. 5 на рисунке 3) for(int pyramidIdK = 0; pyramidIdK < countPyramidsK; pyramidIdK++)

int startPyramidBasePositionK = pyramidIdK * pyramidBaseLengthK; //начальный индекс //верхнего основания пирамиды в сетке getOffsets(pyramidIdK, pyramidBaseLengthK, &currentBaseLengthK, &leftOffsetK,

krightOffsetK);

getBorders(currentBaseLengthK, leftOffsetK, rightOffsetK, kfullBaseLengthK,

&leftPyramidBorderK, &rightPyramidBorderK);

//Копирование полей на видеокарту (п. 6 на //рисунке 3)

create_arrays_on_GFU(Imax,Jmax,fullBaseLengthK); copy_temp_arrays_to_GPU(0,0,leftPyramidBord erK,Imax,Jmax,fullBaseLengthK); copy_constant_to_device();

//Цикл по времени внутри прохода //(п. 7 на рисунке 3)

for (int n = currentTime; n <= durationOfTimePass && n <=nmax; n++) {

//(п. 8 на рисунке 3)

//Ядро для пересчёта компонент магнитного //поля на GPU

kernelH_pyramid1<<<gridSize,blockSize>>>(...); cudaEventSynchronize(syncEvent);

//Ядро для пересчёта компонент электриче-//ского поля на GPU

kenelE_pyramid1<<<gridSize,blockSize>>>(...); cudaEventSynchronize(syncEvent);

//Копирование данных из видеопамяти (п. 9 //на рисунке 3)

copy_arrays_from_GPU_with_offset(currentBas eLengthK, startPyramidBasePositionK, left-OffsetK);

delete_arrays_on_GPU(); currentTime = durationOfTimePass +1;

}

} "'

Сеточные подобласти, соответствующие нижним основаниям пирамид и содержащие начальные данные

для каждого прохода, пересекаются (аналогично показанному на рис. 1, 2 для одномерного и двумерного случаев соответственно), поэтому перед началом очередного прохода необходимо скопировать эти области во временный буфер (реализуется процедурой сге-аге_гетр_/1е1Ж()), чтобы значения не были изменены при копировании из видеопамяти значений, рассчитанных при обработке соседней пирамиды.

3.2. Постановка вычислительных экспериментов

Исследуя условия эффективности применения метода пирамид, проведем две серии вычислительных экспериментов, ограничив объем используемой видеопамяти 128 и 256 Мб, что позволит построить до 16 (в первом случае) и 4 (во втором) пирамид с разной шириной основания и высотой. Целью экспериментов будет определение зависимости ускорения вычислений от данных параметров.

На рис. 4 и 5, в табл. 9 и 10 приведены результаты измерений длительности вычислений.

256x256x256 для 128 Мб

10000

1000

1 2 3 4 5 6 7 8

Рис. 4. Зависимость времени расчёта Тсотр. (с) от высоты (число отсчётов) пирамиды для различной ширины основания (128Мб)

256x256x256 для 256 Мб

ширина 128 ширина 80 ширина 64 ширина 32 ширина 16

2 4 8 16 32 64 Рис. 5. Зависимость времени расчёта Тсотр. (с) от высоты (число отсчётов) пирамиды для различной ширины основания (256Мб)

Как показали результаты экспериментов, минимальное время вычислений достигается при соблюдении баланса использования памяти устройства и числа пересылок. В эксперименте это достигнуто при использовании пирамид с шириной основания 32 и высотой 8 отсчётов при использовании 128 Мб памяти; пирамид с шириной основания 64 и высотой 20 отсчётов при использовании 256 Мб памяти.

Табл. 9. Зависимость времени расчёта (с) от высоты (число отсчётов) пирамиды для различной ширины основания

т г п Объём памяти для расчёта одной пирамиды, Мб Время расчёта, с Ускорение

8 32 1 81 3854 1

8 32 2 84 2091 1,84

8 32 4 91 1253 3,07

8 32 8 105 879 4,38

16 16 1 53 4456 0,86

16 16 2 56 2400 1,61

16 16 4 63 1436 2,68

16 16 8 77 1035 3,72

Табл. 10. Зависимость времени расчёта (с) от высоты (число отсчётов) пирамиды различной ширины основания

т г п Объем памяти для расчёта одной пирамиды, Мб Время расчёта, с Ускорение

2 128 1 246 3403 1

2 128 2 248 1896 1,79

4 80 1 164 3560 0,96

4 80 2 168 1969 1,73

4 80 4 175 1190 2,86

4 80 8 189 820 4,15

4 80 16 216 683 4,98

4 64 1 136 3561 0,96

4 64 2 140 1960 1,74

4 64 4 147 1184 2,87

4 64 8 161 814 4,18

4 64 16 189 673 5,05

4 64 20 202 660 5,16

4 64 32 244 694 4,90

8 32 1 81 3855 0,88

8 32 2 84 2092 1,63

8 32 4 91 1254 2,71

8 32 8 105 880 3,86

8 32 16 133 795 4,28

8 32 20 147 809 4,21

8 32 32 189 956 3,55

8 32 40 216 1062 3,20

8 32 50 251 1177 2,89

16 16 1 53 4457 0,76

16 16 2 56 2400 1,41

16 16 4 63 1436 2,37

16 16 8 77 1035 3,29

16 16 16 105 1054 3,23

16 16 20 119 1121 3,03

16 16 32 161 1464 2,32

16 16 40 189 1682 2,02

При дальнейшем увеличении высоты пирамиды время расчёта начинает увеличиваться за счёт увеличения объема дублирующихся данных, что иллюстрируется и-образной зависимостью времени вычислений от высоты пирамиды. На рис. 4 зависимость имеет линейный характер, однако ограничение памяти в 128 Мб не позволяет более увеличивать высоту пирамиды при указанной ширине основания, в результате чего невозможно достигнуть превышения объема дублирующихся

вычислений над временем пересылки, что хорошо видно при ограничении памяти в 256 Мб.

Заключение

Представленная реализация FDTD с использованием метода пирамид позволяет снизить влияние ограниченного объема памяти графического процессора и использовать преимущество GPU в быстродействии за счёт снижения интенсивности коммуникаций между оперативной и видеопамятью за счёт дублирования арифметических операций.

При реализации вычислений по методу пирамид важно соблюдать баланс между объемом пересылок и задействованным объемом памяти: определяя ширину основания пирамиды при фиксированной высоте, необходимо выбирать максимально допустимую ширину основания, так как это сокращает число пересылок. Расчёты, произведенные с использованием пирамид одинаковой высоты, имеют близкие значения ускорения, что также подчеркивает значительность затрат на пересылку данных по сравнению с затратами на непосредственный расчёт.

Благодарности

Работа выполнена при поддержке гранта РФФИ 14-07-00291-а.

Литература

1. Taflove A. Computational Electrodynamics: The Finite-Difference Time-Domain Method / A. Taflove, S. Hagness. -3th ed. - Boston: Arthech House Publishers, 2005. - 1006 p.

2. Котляр, В.В. Фотонные струи, сформированные квадратными микроступеньками / В.В. Котляр, С.С. Стафе-ев, А.Ю. Фельдман // Компьютерная оптика. - 2014. -Т. 38, № 1. - С. 72-80.

3. Тиранов, А.Д. Коллективное спонтанное излучение в волноводе с близким к нулю показателем преломления /

А.Д. Тиранов, А.А. Калачёв // Известия РАН, серия физическая. - 2014. - Т. 78, № 3. - С. 271-275.

4. Перов, С.Ю. Теоретическая и экспериментальная дозиметрия в оценке биологического действия электромагнитных полей носимых радиостанций. Сообщение 1. Плоские фантомы / С.Ю. Перов, Е.В. Богачёва // Радиационная биология: Радиоэкология. - 2014. - Т. 54, № 1. - С. 57-61.

5. Основы работы с технологией CUDA / А.В. Боресков, А.А. Харламов. - М.: ДМК Пресс, 2010. - 232 с.

6. OpenCL - The open standard for parallel programming of heterogeneous systems / URL: http://www.khronos.org/opencl/.

7. B-CALM - Belgium California Light Machine [Электронный ресурс]. - URL: http://b-calm.sourceforge.net/.

8. FDTD solver [Электронный ресурс]. - URL: http : //www. acceleware.com/ fdtd- solvers.

9. Lamport, L. The parallel execution of DO loops / L. Lamport // Communications of the ACM. - 1974. -Vol. 17(2). - P. 83-93.

10. Вальковский, В.А. Параллельное выполнение циклов. Метод пирамид / В.А. Вальковский // Кибернетика. -1983. - № 5. - С. 51-55.

11. Головашкин, Д. Л. Решение сеточных уравнений на графических вычислительных устройствах. Метод пирамид [Электронный ресурс] / Д.Л. Головашкин, А.В. Кочуров // Современные проблемы прикладной математики и механики: теория, эксперимент и практика (Новосибирск, Россия, 30 мая - 4 июня 2011 г.). - Новосибирск: ИВТ СО РАН, 2011. - URL: http://conf.nsc.ru/files/conferences/niknik-90/fulltext/37858/46076/kochurov_final.pdf.

12. Yee, K.S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media / K.S. Yee // IEEE Transactions on Antennas and Propagation. -1966. - Vol. AP-14. - Р. 302-307.

13. Климов, В.В. Наноплазмоника / В.В. Климов. - М.: Физ-матлит, 2009. - 480 с.

14. Taflove, A. Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Max-wells's equation's / A. Taflove, M. Brodwin // IEEE Transactions of Microwave Theory and Techniques. - 1975. -Vol. mtt-23, Issue 8. - P. 623-630.

Сведения об авторах

Малышева Светлана Александровна, Самарский государственный аэрокосмический университет имени академика С.П. Королева (национальный исследовательский университет), аспирант. Область научных интересов: FDTD-метод, векторные и матричные вычисления, CUDA. Е-mail: [email protected] .

Головашкин Димитрий Львович, доктор физико-математических наук, профессор кафедры прикладной математики Самарского государственного аэрокосмического университета имени академика С.П. Королева (национального исследовательского университета), ведущий научный сотрудник Института систем обработки изображений (ИСОИ) РАН. Область научных интересов: FDTD-метод, векторные и матричные вычисления, математическое моделирование оптических элементов и устройств нанофотоники. E-mail: [email protected] .

Поступила в редакцию 15 февраля 2016 г. Окончательный вариант - 5 апреля 2016 г.

IMPLEMENTATION OF THE FDTD ALGORITHM ON GPU USING A PYRAMID METHOD

S.А. Malysheva 1-2, D.L. Golovashkin 1,2

1 Image Processing Systems Institute, Russian Academy of Sciences, Samara, Russia,

2 Samara State Aerospace University, Samara, Russia

Abstract

In this paper we develop a pyramid method in the context of solving time-dependent Maxwell's equations based on the finite difference time domain (FDTD) approach, which is implemented on a

graphics processing unit (GPU). Application of this method allows the impact of the GPU's limited memory capacity on the computation time to be reduced, which is significant for the FDTD method. Keywords: FDTD, Maxwell's equations, method of pyramids, GPU, CUDA. Citation: Malysheva SA, Golovashkin DL. Implementation of the FDTD algorithm on GPU using a pyramid method. Computer Optics 2016; 40(2): 179-87. DOI: 10.18287/2412-6179-2016-40-2179-187.

Acknowledgements: The work was partially funded by the Russian Foundation of Basic Research Grants No. 14-07-00291-a.

References

[1] Taflove A, Hagness S. Computational Electrodynamics: The Finite-Difference Time-Domain Method. 3th ed. Boston: Arthech House Publishers; 2005.

[2] Kotlyar VV, Stafeev SS, Feldman AYu. Photonic Nanojets Formed By Square Microsteps [In Russian]. Computer Optics 2014; 38(1): 72-80.

[3] Tiranov AD, Kalachev AA. Collective spontaneous emission in a waveguide with close to zero refractive index [In Russian]. Bulletin of the Russian Academy of Sciences 2014; 78(3): 271-275.

[4] Petrov SY, Bogacheva EV. Theoretical and experimental dosimetry in the assessment of biological action of electromagnetic fields portable radio. Message 1. Flat phantoms [In Russian]. Radiation Biology: Radioecology 2014; 54(1): 57-61.

[5] Boreskov AV, Harlamov AA. The basics of working with CUDA technology [In Russian]. Moscow: DMK Press, 2010.

[6] OpenCL - The open standard for parallel programming of heterogeneous systems. Source: (http://www.khronos.org/opencl).

[7] B-CALM - Belgium California Light Machine. Source: (http://b-calm.sourceforge.net/).

[8] FDTD solver. Source: (http://www.acceleware.com/fdtd-solvers).

[9] Lamport L. The parallel execution of DO loops. Communications of the ACM 1974; 17(2): 83-93.

[10] Valkovskii V. Parallel execution cycles. Method of the pyramids [In Russian]. Cybernetics 1983; 5: 51-55.

[11] Golovashkin DL, Kochurov AV. The decision of the grid equations graphical computing devices. Method pyramids. Source: (http://conf.nsc.ru/files/conferences/niknik-90/fulltext/37858/46076/kochurov_final.pdf).

[12] Yee KS. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans Antennas Propag 1966; AP-14: 302-307.

[13] Klimov VV, Nanoplasmonics [In Russian]. Moscow: Fizmatlit; 2009.

[14] Taflove A, Brodwin M. Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Max-wells's equation's. IEEE Transactions of Microwave Theory and Techniques 1975; 23(8): 623-630.

Authors' information

Svetlana Alexandrovna Malysheva, Samara State Aerospace University (National Research University), post-graduate student. Scientific interests: FDTD-method, vector and matrix computations, CUDA. E-mail: s-a-mal@yandex. ru .

Dimitriy Lvovich Golovashkin, Doctor of Physical and Mathematical Sciences, Professor of Samara State Aerospace University (National Research University), leading scientist of Image Processing Systems Institute of RAS. Scientific interests: FDTD-method, vector and matrix computations, mathematical modeling optical elements and devices of nanophotonics. E-mail: [email protected] .

Received February 15, 2016. The final version - April 5, 2016.