CC BY
37
Вычислительные технологии
Том 12, Специальный выпуск 4, 2007
ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ МЕТОДА КОНЕЧНЫХ ОБЪЕМОВ ДЛЯ РЕШЕНИЯ НЕСТАЦИОНАРНЫХ УРАНЕНИЙ МАКСВЕЛЛА НА НЕСТРУКТУРИРОВАННОЙ СЕТКЕ*
Л.Ю. Прокопьева, Ю.И. Шокин, A.C. Лебедев, М.П. Федорук Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: liudmila.prokopeva@gmail.com, shokin@ict.nsc.ru sasa@ict.nsc.ru, mife@ict.nsc.ru
We address parallelization aspects of the finite volume solver for time-dependent Maxwell equations on unstructured grids, based on domain decomposition. The data dependence analysis is first reviewed for derivation of the required graph of dependencies and the scheme of processors' communications. The results of acceleration of computations obtained on the clusters at the HLRS and the ICT SB RAS using MPI for processors' communications are discussed.
Введение
В последнее время в связи с развитием новых технологий в волоконной оптике возросла актуальность численного решения уравнений Максвелла. Создание современных оптоволоконных линий связи порождает широкий спектр задач для научных исследований. Помимо решения задач по оптимизации и проектированию конкретных линий ведутся фундаментальные исследования по использованию новых искусственно созданных материалов с нетривиальными оптическими свойствами. Последние достижения физических экспериментов связаны главным образом с исследованиями сложных композитных конструкций с периодическими включениями, имеющими наноразмеры, В частности, были получены материалы с отрицательным коэффициентом показателя преломления (см., например, работы [1-3] и приведенную в них библиографию). Это, в свою очередь, открывает путь к созданию новых технологий производства современных оптических приборов.
В то же время возможности экспериментальных исследований подобных структур весьма ограничены. В этих условиях особое значение приобретают методы математического моделирования. Однако специфика задач накладывает ограничения и на численные методы. Это прежде всего сложная геометрия, малые (часто порядка длины волны падающего света) размеры включений, а также разрывы решений на разделах сред.
* Работа выполнена при финансовой поддержке Президента РФ (грант № НШ-9886.2006.9) и Российского фонда фундаментальных исследований (грант № 06-01-00210).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
В таких случаях в качестве математической модели не могут быть использованы нелинейные уравнения Шредингера, обычно применяемые для моделирования электромагнитных процессов в световодах, поэтому возникает необходимость в решении полных уравнений Максвелла,
Задачей построения вычислительных алгоритмов для решения нестационарных уравнений Максвелла в областях со сложной формой границ, характерной для этих приложений волоконной оптики, в Институте вычислительных технологий СО РАН начали заниматься еще в 1998 г. Для решения нестационарных уравнений Максвелла на неструктурированных сетках предложен метод конечных объемов (МКО) [4]. Впоследствии был получен алгоритм, позволяющий проводить расчеты для сред с разрывными значениями диэлектрической проницаемости [5], Этот алгоритм, имеющий второй порядок точности, был обобщен для моделирования в средах с дисперсией в рамках моделей Дебая и Друде, Разработанные на данный момент методы на треугольных неструктурированных сетках позволяют адекватно моделировать распространение электромагнитных полей в сложных средах с нетривиальной геометрией границ. На данном этапе проведены одно- и двумерные тесты дифракции плоской монохроматической волны на различных разделах диэлектрик — металл [6], а также на градиентном цилиндре [7], которые хорошо согласуются с аналитическими решениями.
Именно в композитах диэлектрик — металл со сложной геометрией удается получить отрицательное значение коэффициента показателя преломления. Однако для моделирования сколько-нибудь существенного образца такого материала с включениями, имеющими наноразмеры, требуются огромные вычислительные мощности, В одних случаях задачу можно упростить, проводя расчеты на единичном паттерне с постановкой периодических граничных условий. Другой подход состоит в использовании средств параллельного программирования и современных многопроцессорных вычислительных систем,
В этой связи возникла задача распараллеливания разработанного МКО решения уравнений Максвелла для проведения расчетов на современных вычислительных комплексах, Поставленная задача была решена методом декомпозиции расчетной треугольной сетки на вычислительные подобласти для каждого процессора. Обмен необходимых расчетных величин между процессорами реализован с помощью библиотеки MPI, Эффективность кода тестировалась на кластере HLES (High Performance Computing Center, Stuttgart), На малом числе процессоров (до 20), используемом в тестовых расчетах, получено практически линейное ускорение,
В статье описывается опыт распараллеливания МКО на неструктурированной треугольной сетке. Излагается применяемый последовательный алгоритм численного решения уравнений Максвелла и попутно анализируются зависимость по данным и возможность параллелизации, приводится формальная схема полученного параллельного алгоритма и требуемых межпроцессорных обменов. Обсуждаются проблемы декомпозиции области и ее реализации с использованием теневых (или фиктивных) ячеек (halo cells, ghost cells, shadow cells), Приводятся расчеты ускорения параллельного кода, реализующего полученный алгоритм, написанный на языке С++ с использованием интерфейса MPI,
Опыт проделанной работы и полученные положительные результаты могут оказаться полезными для начала работы над основанными на декомпозиции области параллельными версиями других численных алгоритмов и приложений.
1. Параллельная версия МКО для решения уравнений Максвелла
Распространение электромагнитных волн в немагнитной оптической среде (р = 1, ^уВ = 0) описывается уравнениями Максвелла, которые запишем в следующем виде (см, например, [8]):
ЭР
~дГ ав
!к
- гоШ
0, Б(ш) = £г*НЕ
+ го1Е &уБ
0, 0,
В = Н,
а1уВ = о.
(1)
Здесь дисперсионные свойства среды описываются зависимостью £*(ш>), которая для среды с относительной диэлектрической проницаемостью £г и конечной проводимостью а, описываемой в классической системе уравнений Максвелла, имеет вид
£*(о>) = £г + а/гш£0. Важными для приложений являются модели Дебая и Друде:
£*(ш) = +
Х1
а
1 —
гш£0
(2)
(3)
£*(ш) = —
и(и + гГ)'
(4)
где — высокочастотная диэлектрическая проницаемость; Х1 — линейная восприимчивость среды; а — проводимость среды; шр — циклическая плазменная частота; Г — коэффициент затухания.
Физически модель Друде применяется для металлических сред, что соответствует отсутствию первого слагаемого в разложении по 1/и. С математической же точки зрения модель Друде — лишь частный случай модели Дебая, поэтому конечную для численного решения систему уравнений рассмотрим лишь для модели Дебая, Используя Б(ш) = £* (ш)Е и выражение (3), перейдем во временную область с помощью обратного преобразования Фурье, В результате получим
БШ = £ооЕ(*) + — / е-(*-т)/*°Е(т)с1т + — I Е(т)с1т.
¿о У £о У
пп
(5)
После переобозначения переменной В ^ £юЕи подстановки ее в уравнение (1) получим окончательно систему уравнений, описывающую распространение электромагнитных полей в дисперсионных средах:
д Б
— гоШ = —д,
— + го!Е = 0, дЬ
&Уб а1уВ
0, 0,
Б
В
Н,
£
£
где
1-П
е-(*-г)/*о+ ^ + ^ I Е(т)с1т.
(7)
0 0 Полученная система уравнений теперь учитывает дисперсионные свойства материалов и адекватно описывает электромагнитные явления для широкого круга современных материалов, применяемых в задачах волоконной оптики.
Далее будем использовать матричную запись уравнений для двумерного случая (осесимметричный случай, ТМ-, ТЕ-волны) в дивергентной и недивергентной формах соответственно:
д и д д
(8)
ди д , тт д — + —л и + —А2 и = -И;
дг дх1 дх2
дУ дУ дУ
-Ц.
(9)
Например, для распространения ТМ-волны в металлической среде, описываемой моделью Друде, матрицы уравнений (8) и (9) имеют вид
¿1
и
В1
0 0
1
0 -1 \ 00
00
¿2
/
( 0 1 0 \
.1 о о V 0 0 0 /
Вз Н1 Н2
и
О г ^
0 0
/ 0 0 1 \ 0 1 0
0 0 0 , В2 = 1 0 0
V -1 0 0 1 0 0 0
(10)
У
Ез Н1 Н2
Ц
\
0
0 0
/
(Н)
Для построения численной схемы метода конечных объемов в расчетной двумерной области строится неструктурированная треугольная сетка, учитывающая геометрию внешних границ области и раздела сред (см., например, рис, 1),
Уравнение (8) интегрируется по каждой треугольной ячейке Д» и по формуле Гаусса — Остроградского преобразуется к контурному интегралу
J \Jdxidx2 + ^ J АШГ = - J И(1х
Дг к=1 гк Д,
где Гк — ребра ячейки; п = (п1, п2) — внешняя нормаль; А = п1А1 + п2А2,
4
í
í
Рис. 1. Сетка в области с вырезом, 2481 ячейка
Определим значения поля ИП в барицен тре хь треугольного эле мента Дг в момент времени Ьп = ит, где т — шаг по времени. Кроме того, обозначим приближенное значение потока Аи в центре &-го ребра г-го треугольн ика Ег,к, длину ре бра /г>к, площадь г-го элемента Бг, Аппроксимируя уравнение (12) в момент времени Ьп = ит в центре г-й ячейки, получим явную схему
т тп+1 т тп 3
т
к=1
Б^-И".
(13)
На каждом временном слое вычисления выполняются явно. Обусловлено это тем, что потоки зависят от данных в соседних ячейках, что обусловливает основные трудности при распараллеливании.
Поток Ег,к предлагается вычислять как решение одномерной относительно нормали к &-му ребру г-го треугольника задачи Римана, записанной в дивергентных переменных:
Е-
А+И + А-И
Я 5
(14)
где
А = Б-1diag
-(Ак±|Ак|)^,
Б
А; Ак — собственные числа матрицы А; и^,Я — значения поля и, аппроксимированные на ребро из левого и правого примыкающих треугольников.
Например, для рассмотренного случая ТМ-волны в модели Друде матрицы выглядят следующим образом:
А+
/ 1
П2 -П1
\/е п|
1 п2 щп2
2 £ л/е л/е
пг П\П2 и\
V £ Те
\
А-
/
/ 1
и -и1
\/е _ п|
1 п2 щп2
2 £ л/е у/ё
пг П\П2 и1
V £
\
2£ь£Д £ь + £я .
А
раллельпой версии алгоритма матрицы А± для каждого ребра могут быть вычислены только один раз и хранятся локально в каждой подобласти.
Аппроксимация же значений поля из центра треугольников на ребро, т, е, вычисление Ид и Иь, требует дальнейших пояснений (рис, 2), Для аппроксимации будем использовать недивергентные переменные Уд и V ь (Ид,ь = ТУд,ь, например, для нашего случая Т = diag(£, 1,1)). Разложив поле V по формуле Тейлора в окрестности точки (¿га,жЬ), получим
у,-/*= V к, д+Е к - 4.) + (15)
%=1 %
или, после подстановки производной по времени из уравнения (9),
%=1 %
Таким образом, для вычисления потока через каждое граничное ребро потребуется явно пересылать значение Уь(УД), аппроксимированное из внешнего смежного треугольника, принадлежащего соседней области, С этой целью внутри каждой подобласти для каждого внутреннего треугольника эти значения должны быть вычислены для каждого из трех ребер. Согласно (16) при этой аппроксимации требуется знать градиенты дИ/дх% поля в каждом внутреннем треугольнике области. Здесь, как мы увидим далее, потребуются дополнительные межпроцессорные обмены.
Перейдем к вычислению градиентов поля в каждой подобласти. Эта процедура осуществляется в несколько этапов и является ключевым местом алгоритма при работе с разрывными профилями диэлектрической проницаемости е. В вычислении градиентов помимо значений поля в непосредственно примыкающих треугольниках участвуют и значения поля в следующих за соседними треугольниках. Это означает, что, вообще говоря, для параллельной реализации такого вычисления потребуется пересылать значения поля не только в элементах, непосредственно примыкающих к границам областей, но также и следующих за ними, В данном случае этого удается избежать за счет дополнительных обменов после каждого этапа вычисления градиентов.
Этап 1. Пусть поле V известно в центре каждого треугольника Д%, Вычислим предварительные градиенты поля в каждом треугольнике, используя значения поля в трех
Рис. 2. Две соседние ячейки Дь и Дд треугольной сетки с барицентрами хЬ и хД и общим ребром, середина которого находится в точке х^
Рис. 3. Обмен значениями поля V для вычис- Рис. 4. Обмен значений, аппроксимируемых ления предварительных градиентов в вершины из внешних треугольников. Пере-
сылать достаточно сумму для последующего усреднения
Процессор О Процессор М-1
Рис. 5. Блок-схема параллельного МКО решения уравнений Максвелла
соседних треугольниках (случай, когда ячейка имеет менее трех соседей, здесь опускается), Эта процедура может быть осуществлена параллельно, если для каждого граничного ребра значение поля во внешнем фиктивном треугольнике будет передано заранее (рис, 3),
Эт,ап 2. С помощью найденных предварительных градиентов аппроксимируем значения поля из центров ячеек в вершины. Эта процедура может быть выполнена локально в каждой области, поскольку для этого требуется лишь значение поля и градиентов в каждом внутреннем треугольнике области. Однако итоговое значение поля в вершине берется как среднее из значений, полученных аппроксимацией из каждого инцидентного треугольника. Таким образом, передавать из каждой области нужно только одно суммарное значение поля, при этом предполагается что для каждой граничной вершины глобальное число инцидентных треугольников известно (рис, 4),
Эт,ап 3. Внутри каждой области локально находим окончательные градиенты поля в центре каждого треугольника по значениям поля V в трех вершинах, вычисленных на предыдущем этапе.
Таким образом, алгоритмически исследуемый МКО для решения уравнений Максвелла можно считать распараллеленным. Блок-схема параллельной версии алгоритма приведена на рис. 5. Алгоритм требует трех сеансов обменов на каждом шаге по времени, Технические подробности в данной статье опущены, некоторые проблемы реализации обсуждены.
2. Декомпозиция области и тесты ускорения параллельного кода
Для проведения тестовых расчетов выбиралась квадратная расчетная область с числом треугольников на сторонах 80, 160, 320, 640. Тестовые расчеты выполнялись для задачи распространения электромагнитной ТМ-волны в диэлектрическом слое [9]. По про-
Рис. 6. Пример декомпозиции области горизонтальными полосами. Число треугольников на стороне квадрата 20
Рис. 7. Другой пример декомпозиции той же области
цессорам ячейки делились геометрически — горизонтальными полосками и в форме "паутины" (рис, 6 и 7),
Сетка каждой из подобластей хранилась локально на некотором процессоре. Ключевым фактором при организации и программной реализации обменов являются теневые ячейки (halo cells), которые позволяют сохранить все процедуры последовательной программы, Для этого к каждой сетке надо достроить слой теневых ячеек и обеспечить своевременную передачу в них значений, необходимых для каждого текущего расчета. Дополнительные технические трудности возникают лишь с синхронизацией фиктивных и соответствующих им внутренних ячеек некоторой области.
Расчет ускорения производился на кластере HLRS, технические характеристики оборудования приведены в таблице. Для сетки с достаточно большой размерностью при горизонтальном делении получено практически линейное ускорение. Проблемы со вторым методом деления связаны с балансировкой ячеек по процессорам (рис, 8 и 9), Там, где балансировка близка к идеальной (для 17 процессоров), ускорение по-прежнему остается линейным. Таким образом, на малом числе процессоров, вплоть до 20, параллельная версия рассматриваемого МКО решения уравнений Максвелла обеспечивает линейное ускорение при условии балансировки и достаточно большой размерности. Для малой размерности (80x80 треугольников) плохое ускорение, по-видимому, связано с существенной долей обменов в вычислениях.
Характеристика используемого оборудования
HLRS II BT CO PAH
Number of nodes 210 nodes 4 nodes
A node Dual 3.2 GHz Intel Xeon EM64T Dual 3.06 GHz Intel Xeon
Memory 1Gb to 8Gb Memory 2Gb Memory
Network Infiniband Gigabit Ethernet
MPI Voltaire Grid Director ISR 9288 MPICH
Voltaire MPI
CC Intel Compiler 9.0 (EM64T) Intel Compiler (8.1, 9.1)
QMS Open PBS Sun Grid Engine 6.0 nlge
Рис. 8. Расчет ускорения для горизонтальной декомпозиции
Рис. 9. Расчет ускорения для декомпозиции, представленной на рис. 7
11 10
9
t 8 I-T 7
II
6
3 Т5
0) 5
0) ^
(Л 4
3
Y080
Lin
Hor_080_HLRS Web_080_HLRS Hor_080_ICT Web 080 ICT
J_I_I_I_I_I_I_I
23456789 10
Number of Processors
ii 10
9
tf 8 I-" 7 II
6
3 T5
0) 5
0)
Ci
(Л 4
3
Y160
- Lin
- Hor_160_HLRS
- Web_160_HLRS
-a- Hor_160_ICT
-0- Web_160_ICT
23456789
Number of Processors
2
Рис. 10. Ускорение на кластерах HLRS и ИВТ СО РАН. С увеличением размерности ускорение выходит на линейный рост на обоих машинах
Интересно сравнить результаты ускорения кода на кластере HLRS с ускорением на менее мощном, но архитектурно схожем кластере ИВТ СО РАН (см. таблицу). Из рис. 10 видно, что с ростом размерности задачи ускорение становится менее чувствительным к разнице в пропускной способности сети (Infiniband и Gigabit Ethernet) и выходит на линейный рост для обоих вычислительных кластеров.
Заключение
Рассмотрены алгоритмические аспекты распараллеливания метода конечных объемов на неструктурированных сетках. Полученный код с использованием интерфейса MPI демонстрирует практически линейное ускорение на кластерных архитектурах (см. таблицу) для достаточно большой размерности и рассматриваемого числа процессоров (Np < 20).
Дальнейшие направления работы над параллельным кодом связаны с автоматизацией декомпозиции сетки, применением алгоритмов на графах вместо геометрического
деления, используемого на предварительном тестовом этапе.
Авторы выражают благодарность Центру высокопроизводительных вычислений
III.HS (Germany, Stuttgart) за предоставленные для вычислений ресурсы и лично
М, Рэшу и Т. Бёнишу,
Список литературы
fl] Smith D.R., Pendry J.B, Witshire M.C.K. Metamaterials and negative refractive index // Science. 2004. Vol. 305, N 3, Aug. P. 788-792.
[2] markos P, Soukoulis C.m. Trasnsmission properties and effective electromagnetic parameters of double negative metamaterials // Optics Express. 2003. Vol. 11, N 7, April. P. 649-661.
[3] Kildishev A.V. et al. Negative refractive index in optics of metal-dielectric composites // J. Opt. Soc. Am. B. 2006. Vol. 23, N 3, March.
[4] Fedoruk M, Münz C.-D, Omnes P, Schneider R. A Maxwell-Lorentz Solver for Self-Consistent Particle-Field Simulations on Unstructured Grids. Karlsruhe: Forschungszentrm Karlsruhe GmbH, 1998.
[5] Лебедев A.C., Федорук М.П, Штырина O.B. Конечно-объемный алгоритм решения нестационарных уравнений Максвелла на неструктурированной сетке // Журн. вычисл. математики и мат. физики. 2006. Т. 46, № 7. С. 1302-1317.
[6] Kildishev A.V, Chettiar U. Realistic Models of Metal Nano-Scaterers at Optical Frequencies. Interim Technical Report. School of ECE, Purdue Univ., June 2004.
[7] Котляр B.B, Личманов M.A. Дифракция плоской электромагнитной волны на градиентном диэлектрическом цилиндре // Числ. методы компьютерной оптики. 2003. Т. 25. С. 11-15.
[8] Sullivan D.M. Electromagnetic Simulation Using the FDTD Method. N.Y.: The Institute of Electrical and Electronics Engineers, Inc. 2000.
[9] Advances in Computational Electrodynamics. The Finite-Difference Time-Domain Method / A. Taflove (ed.). Boston: Artech House, 1998.
Поступила в редакцию 31 июля 2007 г.
