Разностное решение волнового уравнения на графических процессорах с повторным использованием попарных сумм дифференциального шаблона

Разностное решение уравнений Максвелла (FDTD-метод [1]) на протяжении последних 15 лет является незаменимым инструментом математического моделирования в быстроразвивающихся областях естествознания, изучающих генерацию, распространение и приём электромагнитного излучения: нанофотоники [2], наноплазмоники [3], теории фотонных кристаллов [4] и др. Характеризуясь простотой программной реализации, общностью математической модели, естественностью подхода к её решению, FDTD-метод вместе с тем известен высокой требовательностью к системным ресурсам вычислительного устройства, препятствующей его широкому распространению за пределами сообщества пользователей суперкомпьютеров.

Решение указанной проблемы традиционно ищут в упрощении математической модели либо в переходе к вычислениям на распространённых графических процессорах (GPU). В качестве примера первого подхода уместно сослаться на разностное решение уравнения Даламбера [1, 5, 6], применение которого при известных условиях [7] (TM-случай) позволяет получать решение без потери точности по сравнению с классическим FDTD-методом, при снижении на треть требований к объёму памяти ЭВМ. Следование второму подходу [8, 9, 10] позволяет задействовать богатые вычислительные ресурсы современных видеокарт, содержащих тысячи исполнительных устройств (CUDA-ядер) и находящихся в распоряжении большинства исследователей. Признанным недостатком GPU является небольшой объём видеопамяти по сравнению с кластерными решениями [9], что делает актуальным сочетание обоих упомянутых подходов, а именно – разностное решение уравнения Даламбера на графических процессорах.

Признавая несомненный успех применения GPU при разностном решении уравнений Максвелла (так, в [8] иллюстрируется, как использование пакета B-CALM, ориентированного на видеопроцессор, позволяет в 30 раз сократить длительность расчётов по сравнению с пакетом MEEP, использующим обычный кластер), необходимо отметить определённое замедление в развитии данного направления вычислительной математики. Переход к современной методике построения векторных алгоритмов, основанной на использовании «длинных» векторов, позволил лишь незначительно (в 1,4 раза [11]) ускорить вычисления по сравнению с B-CALM. Это обстоятельство связано с особенностями дифференциального шаблона классической схемы Yee [1], ограничивающими выбор варианта векторизации вычислений алгоритмом на основе операции gaxpy (general x plus y [12]), уступающим по ускорению алгоритму на основе повторного использования попарных сумм [11]. Вместе с тем, как будет показано далее, использование дифференциального шаблона для уравнения Даламбера [13] освобождает от упомянутого ограничения, что также свидетельствует об актуальности выбранного направления исследований.

• 1.    Векторный алгоритм на основе повторного использования попарных сумм

Заботясь в первую очередь о демонстрации идеи построения векторного алгоритма и сравнении с FDTD-методом, авторы выбрали для их иллюстрации классическое двумерное однородное линейное нестационарное волновое уравнение [7].

д 2U _ 2 ( д 2U d 2U ) =c + д t2      ( д y2 д z2 J,

где U ( t , y ,0) – напряженность электрического поля, c – скорость света в среде (в однородном пространстве – константа, в неоднородном – функция координат). На краях вычислительной области D (0 ≤ t ≤ T , 0 ≤ y ≤ L y , 0 ≤ z ≤ L z ) устанавливаются граничные условия Дирихле U ( t , y ,0) =0 и U ( t , y , L z ) =0, U ( t , 0, z ) =0 и U ( t , L y , z ) =0, соответствующие электрической стенке. В качестве начального условия принимается отсутствие поля в области при t ≤ 0 ( U (0, y ,z) и ∂ U (0, y , z) /∂ t =0). Далее, при t >0 поле появляется из «жёсткого» источника ( hard source , [1]), расположенного в произвольном месте D .

Для конструирования сеточной области воспользуемся известной областью Yee для FDTD-метода [1], исключив узлы, в которых определялись сеточные аналоги проекций магнитного поля Hy и Hz. Получим Dh как {(tn, ym, zk): tn = nht, n = 0, 1, .., K = T/ht, ym = mhy, m = 0, .., M = Ly / hy, zk = khz, k = 0, .., N = Lz / hz}.

Следуя за [13], заменим в (1) частные производные их разностными аналогами и запишем для квадратной области ( L y = L z , M = N , h y = h z = h, h t = τ) следующую хорошо известную явную разностную схему:

U mi - 2 u^u^

τ

-

c ²

n

U m -1, k

-

n

m,k

n

+U m+ 1, k

h ²

+

n

U m,k -1

-

2 U

n m,k

h ²

n

+U m,k+ 1

которую для производства расчётов удобно представить в итерационном виде:

n+1           n          n            n        nn m,k           m,k-1      m-1,k        m,k      m+1,km,k+

+ 2 Unk — U n_m -1 , a = c-^.

m,k m,k2

h

На той же сеточной области соответствующим образом дискретизируются краевые, начальное условия и «жёсткий» источник.

При конструировании длинновекторного алгоритма решения (2) выберем стратегию повторного использования попарных сумм, как принесшую успех в [11] при решении уравнения теплопроводности. Если для схемы Yee в силу особенностей строения сеточной области это невозможно, то после перехода к уравнению Даламбера, – вполне. Ведь дифференциальный шаблон на n -м слое по времени для (2) основан на паттерне «крест», аналогичном для разностного уравнения теплопроводности (рис. 1 из [11]), в то время как по сравнению со схемой Yee это совершенно разные шаблоны. Интересно при этом отметить, забегая вперёд, что решение (2) совпадет с результатами вычислений по схеме Yee в любом знаке и совершенно отличается от решения уравнения теплопроводности в силу разности математических моделей.

Разумеется, уравнение (2) можно переписать и в канонической форме А.А. Самарского, организовав векторизацию вычислений по ней, как это было проделано в [11] для уравнений теплопроводности и Максвелла, однако там же была показана предпочтительность декларируемого подхода с использованием попарных сумм на примере уравнения теплопроводности.

Пользуясь нотацией, введённой в [11,12], запишем векторный алгоритм решения (2) в следующем виде:

Алгоритм

% Заполнение вектора T попарными суммами:

• 1.    T (2: N ( N -1))= V(2:N ( N -1))+ V ( N +1: N*N -1);

• 2.    Z = V ;

• 3.    V(N +2: N(N -1)-1) = a( T (2: N(N -2)-1)+

• 4.    V1 = Z ;

% Сохранение значений Un

% Подсчёт значений для следующего временного слоя U n+1_:

+ T ( N +3: N ( N -2)+ N +1))+a 1 V (( N +2: N ( N -1)-1)+

+ V1 ( N +2: N ( N -1)-1); a 1 =2-4a

%. Восстановление граничных значений.

Вектор nn   n  nn   n   n  n   n v \у 0,0 , 0,1,..,   0,N ,   1,0 , 1,1,..,   1,N ,..,   N,0 ,   N,1,..,   N,N )

формируется построчным «разворачиванием» значений сеточного аналога функции U на временном слое n , вектор V 1 – на слое n– 1. Вектор Z используется для сохранения значений сеточной функции при переходе на следующий слой, когда решение на слое n+ 1 (третья строка Алгоритма) пишется поверх решения на слое n .

Ускорение вычислений по данному алгоритму по сравнению с «коротковекторным» подходом из [14], предложенным для организации вычислений по FDTD-методу, должно обеспечиваться за счёт двух факторов.

• 1.    Перехода к «длинным» векторам, что позволит более полно загрузить вычислительные мощности видеопроцессора. Действительно, если следовать [14] и векторизовать вычисления вдоль одного выбранного пространственного направления сеточной области, то придётся работать с векторами длины ~ N , в предложенном алгоритме производятся операции над векторами из ~ N ² значений.

• 2.    Повторного использования попарных сумм дифференциального шаблона, записываемых в вектор T на первом шаге алгоритма. Например, при вычислении U ц¹ по (2) используется U "2 + U 22,1 , эта же сумма необходима при нахождении значения U " + 21 .

2. Вычислительные эксперименты

При его отыскании нет нужды пересчитывать её второй раз, как это принято в [14].

Проблема восстановления граничных значений (строка 4) и её программное решение аналогичны для разностного решения уравнений теплопроводности и Даламбера. Её обсуждение для первого уравнения приводится в [13].

В отличие от алгоритма для решения уравнений Максвелла из [11], в предложенном снижено количество арифметических операций с векторами вдвое, что позволяет надеяться на двукратное ускорение вычислений. Вместо 4 операций saxpy [12], используется одна, вместо 4 операций сложения – 3 такие операции. При этом производится дополнительное копирование векторов (строки 2 и 5 алгоритма) и восстановление граничных значений, что может снизить выигрыш во времени вычислений.

Предлагаемая работа является завершающей в цикле статей, начатом публикацией [11] и продолженном в [15]. В силу чего программно-аппаратная база для постановки вычислительных экспериментов (видеокарта NVIDIA GeForce GTX 660Ti, центральный процессор Intel Core i7-3770, операционная система Debian Wheezy с установленными драйверами CUDA Toolkit 4 и компилятором gcc), приёмы программной реализации (использование библиотеки CUBLAS для исполнения векторных операций) и па-

раметры вычислительных экспериментов остаются прежними.

На рис. 1 и в табл. 1 приведены сравнения длительности вычислений по различным реализациям разностного решения уравнений Максвелла и Даламбера. Каждое значение получено в результате усреднения по 250 замерам.

Рис. 1. Сравнение длительности вычислений по B-CALM [8] (сплошная линия) и длинновекторными алгоритмами: из [11] (пунктирная линия), основанным на векторизации операции gaxpy, и предложенным в настоящей работе (точечная линия), основанным на повторном использовании попарных сумм дифференциального шаблона

Табл. 1. Сравнение длительности вычислений по B-CALM (T b ), алгоритму из [11] (T 1 ) и предложенному алгоритму (T 2 ). Ускорения S 1 = T 1 / T 2 , S 2 = T b / T 2

N×N	T 1 , с	T 2 , с	S 1	T b , с	S 2
50×50	2,53	2,47	1,10	2,41	0,96
100×100	6,01	2,71	2,21	4,36	1,6
150×150	6,92	2,80	2,47	6,87	2,45
200×200	7,87	3,03	2,59	7,92	2,61
250×250	8,73	3,36	2,58	11,53	3,43
300×300	9,51	3,67	2,6	13,31	3,62
350×350	10,27	3,95	2,61	14,37	3,64
400×400	11,20	4,31	2,59	15,68	3,63
450×450	12,12	4,66	2,60	16,30	3,49
500×500	13,83	5,90	2,34	17,45	2,96
650×650	18,03	8,09	2,22	22,64	2,8
800×800	20,42	10,3	2,10	24,14	2,34
1000×1000	23,86	13,05	1,82	28,50	2,18

Как и ожидалось, длительность вычислений по сравнению с алгоритмом из [11] снижена. Незначительно для небольшой сеточной области ( S 1 = 1,1 при N = 50), существенно для области средних размеров ( S 1 = 2,61 при N = 350) и в соответствии с прогнозом (двукратное уменьшение количества векторных операций) для крупных областей ( S 1 =2,1 и 1,82 при N =800 и 1000).

Для средних и крупных областей более существенным оказался выигрыш авторского алгоритма по сравнению с пакетом B-CALM (табл. 1). К сожалению, программный код пакета, разработанного в Калифорнийском технологическом институте, явля- ется достаточно запутанным, что затрудняет поиск причин достигнутого выигрыша.

Заключение

В настоящей работе предлагается векторный алгоритм разностного решения уравнения Даламбера, предназначенный для моделирования дифракции электромагнитного излучения в двумерном случае. Данным алгоритмом дополняется материал авторских публикаций [11, 15] и завершаются исследования по разработке векторных алгоритмов решения сеточных уравнений различных типов явных разностных схем (двух- и трёхслойных) для различных дифференциальных уравнений (переноса, теплопроводности и волнового), позволяющих ускорить вычисления на GPU по сравнению с популярными пакетами OpenCurrent (компании NVidia) и B-CALM.

Касаясь рекомендаций по использованию алгоритмов из упомянутого цикла работ, авторы в первую очередь желают обратить внимание читателя на решение обратных задач дифракционной нанофотоники [2]. Действительно, использование для этой цели стохастических вычислительных процедур (например, генетического алгоритма из [2]) потребует многократного (десятки тысяч и более) решения прямой задачи дифракции. Тогда трёхкратное ускорение вычислений по сравнению с пакетом B-CALM позволит получить выигрыш по времени не в несколько секунд, как в табл. 1, а значительно более существенный. К тому же в ходе многочисленных вычислительных экспериментов была установлена (и подтверждена разработчиками B-CALM в личной переписке) невозможность использования указанного пакета для решения обратных задач в силу его невысокой надежности.

Работа выполнена при частичной поддержке грантов РФФИ 14-07-00291-а и 16-47-630560-p_a.