автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Разработка комплекса программ решения электродинамических задач с использованием массивно-параллельных вычислительных систем

Московский государственный университет им. М.В. Ломоносова. Факультет вычислительной математики и кибернетики

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

систем

05.13.18 - Математическое моделирование, численные методы и комплексы

На правах рукописи

Семенов Алексей Николаевич

программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

5 ДЕК 2013

Москва - 2013

005542030

005542030

Работа выполнена на кафедре автоматизации научных исследований факультета вычислительной математики и кибернетики Московского государственного университета имени М.В. Ломоносова.

Научный руководитель: к.ф.-м.н., доцент, кафедры автоматизации

научных исследований факультета ВМК МГУ имени М.В. Ломоносова, Смирнов Александр Павлович Официальные оппоненты: д.ф.-м.н., профессор кафедры математической физики факультета ВМК МГУ имени М.В. Ломоносова,

Ильинский Анатолий Серафимович; д.т.н., профессор кафедры 406 факультета РЛА МАИ,

Гринев Александр Юрьевич Ведущая организация: Московский государственный технический

университет радиотехники, электроники и автоматики

Защита состоится 25 декабря 2013 г. в 15 ч. 30 мин. на заседании диссертационного совета Д 501.001.43 при Московском государственном университете имени М.В. Ломоносова, расположенном по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет вычислительной математики и кибернетики, ауд. 685.

С диссертацией можно ознакомиться в научной библиотеке МГУ имени М.В.Ломоносова.

Автореферат разослан

2013 г.

Ученый секретарь диссертационного совета,

доктор физико-математических наук, профессор ^п Л^п ^ Захаров Е.В.

Общая характеристика работы

Актуальность работы. При моделировании распространения электромагнитных волн в радиотехнике [10],[11],[5], микроэлектронике [15],[14], задачах нанооптики [16],[17], [7] и биодиагностики [12], [13], требуется решение полной трёхмерной системы уравнений Максвелла на сетках больших размеров. При реализации на вычислительной системе размер таких сеток составляет порядка терабайта и более. Большие сетки необходимы для описания неоднородностей, характерный размер которых значительно меньше длины волны, а общей размер исследуемой системы имеет размер равный десяткам длин волн. Для решения таких задач необходимы вычислительные системы с большим объёмом памяти и высокой производительностью: суперкомпьютеры с массивно-параллельной архитектурой. Для их использования требуются специализированные комплексы программ, способные эффективно использовать терабайтиые объёмы оперативной памяти и более чем тысячи процессоров вычислительной системы. При реализации комплекса программ необходимо также учитывать специфику архитектуры конкретной вычислительной системы. Комплекс должен быть многофункциональным и рассчитанным на решение задач как в ограниченной, так и в не ограниченной области.

Определение электромагнитных полей в произвольной, неоднородной диэлектрической среде - это важный практический предмет исследования волновых эффектов, в том числе и на микро и нано уровне. Аналитические решения получены, как правило, для объектов простой формы. Для объектов сложной структуры необходимо прибегнуть к численному моделированию задачи. Численное решение полной системы уравнений Максвелла позволяет получить значения компонент электромагнитного поля для последующего исследования и анализа электромагнитных свойств моделируемых структур.

Одним из наиболее распространённых способов численного решения пол-

ной системы уравнений Максвелла является явный метод конечных разностей во временной области, в иностранной литературе finite-difference timedomain (FDTD) [1]. Метод основан на использовании сетки Йи [3], в узлах которой располагаются компоненты электромагнитного поля. Основываясь на уравнениях Максвелла в интегральной форме, можно построить их конечно-разностную аппроксимацию. Поскольку FDTD-метод решает задачу во временной области, он позволяет решать импульсные и негармонические задачи и получить результат для широкого спектра длин волн за один расчёт. Это необходимо при решении задач, в которых неизвестны резонансные частоты, или в случае моделирования широкополосных сигналов. Также FDTD-метод позволяет проследить временную эволюцию распространения волны в моделируемом объёме, что может быть необходимо для детального исследования процесса формирования устойчивой картины распределения электромагнитных волн. Метод позволяет непосредственно моделировать краевые эффекты и эффекты экранирования и удобен при задании анизотропных, дисперсных и нелинейных сред.

При численном решении FDTD-методом задач в неограниченном пространстве на границах расчётной области чаще всего используется искусственный поглощающий слой PML (Perfectly matched layer) [8], [6]. Рассеянная волна, уходя на бесконечность, поглощается анизотропным слоем на границе и, при правильном выборе параметров поглощающего слоя и разностной сетки, отражение от PML слоя минимально [4],[2].

Величина шага дискретизации по пространству должна быть значительно меньше исследуемых длин волн и размера исследуемой структуры. Это может потребовать сеток с маленьким шагом по пространству, что означает большие затраты памяти и возросшее время расчёта. При моделировании дисперсных и анизатропных сред также растут необходимые вычислительные ресурсы. Актуальной становится задача создания программного комплекса

4

на высокопроизводительных массивно-параллельных вычислительных системах с большим объёмом оперативной памяти, который позволит эффективно решать класс электродинамичских задач, описываемых полной системой уравнений Максвелла.

Существует ряд комплексов, реализующих параллельный FDTD-метод на кластерных вычислительных системах([18], [19]). Однако, использование таких комплексов, как правило, рассчитано на массивно-параллельные вычислительные системы общей архитектуры, что делает невозможным эффективно использовать их на специализированных суперкомпьютерах и полностью задействовать их вычислительные мощности. На факультете ВМК Московского государственного университета имени М.В. Ломоносова установлен суперкомпьютер IBM BlueGene/'P, имеющий 8192 процессоров и 4 терабайта оперативной памяти. Разработанный в диссертации программный комплекс был эффективно реализован для этого суперкомпьютера и вычислительных систем схожей архитектуры.

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

Научная новизна. Разработан программный алгоритм для решения крупномасштабных электродинамических задач, который с помощью гибридной технологии MPI и ОрепМР может эффективно выполнятся на современных массивно-параллельных вычислительных системах. На основе предложенного алгоритма создан программный комплекс на многоядерных вычислительных системах и суперкомпьютерах серии IBM Blue Gene/P. Показано, что комплекс имеет хорошую масштабируемость и может эффективно

5

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

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

На защиту выносятся следующие основные результаты и положения:

• Разработан программный алгоритм для решения электродинамических задач, использующий гибридные технологии MPI и ОрепМР. Создан многофункциональный комплекс программ на многоядерных вычислительных системах и суперкомпьютере IBM BlueGene/P, позволяющий решать широкий круг электродинамических задач. Показано, что созданный комплекс имеет хорошую масштабируемость и может эффективно выполняться на массивно-параллельных системах, оперируя с терабайтами вычислительных данных.

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

• С помощью комплекса впервые решена задача вычисления добротности

открытого резонатора Фабри-Перо со вставкой из метаматериала, ранее проводившегося только для не субволновых размеров резонатора.

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

Апробация работы. Основные результаты диссертации докладывались на:

• Международной конференции "Progress in Electromagnetics Research Symposium"(Stockholm, Sweden, 2013 год)

• Международной конференции "Progress in Electromagnetics Research Symposium" (Moscow, Russia, 2012 год)

• Международном научном семинаре "Inverse Problems and Large-Scale Modeling with Applications in Electromagnetics"(Karlstad, Sweden, 2012 год)

• VI Всероссийской конференция "Радиолокация и радиосвязь" Российская академия наук Институт радиотехники и электроники им. В. А. Котелышкова (Москва, 2012 год)

• Объединенном Фельдовском семинаре по электродинамике и антеннам. (Москва, 2012 год) Российская академия наук Институт радиотехники и электроники им. В. А. Котелышкова

• XIII Всероссийской школы-семинара "Физика и применение микроволн"(Зв МГУ им. М.В. Ломоносова, 2011 год)

• Международной суперкомпьютерной конференции "Научный сервис в сети Интернет: экзафлопсное будущее"(Новороссийск, 2011 год)

• Научной конференции "Ломоносовские чтения"(Москва, МГУ им. М.В. Ломоносова, 2011 год)

Публикации. Материалы диссертации опубликованы в 9 печатных работах, из них 6 статей в журналах списка ВАК ([1]-[6]) и 3 тезиса докладов ([7]-[9]).

Структура и объем диссертации Диссертация состоит из введения, 3 глав, заключения и библиографии. Общий объем диссертации 100 страниц, из них 89 страниц текста, включая 23 рисунка. Библиография включает 60 наименований на 7 страницах.

Содержание работы

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

В первой главе диссертации описываются численные методы, используемые при создании комплекса программ численного моделирования уравнений Максвелла ЕМ\УЗо1уегЗО. Описываются разработанные параллельные алгоритмы численного решения, используемые в комплексе. Представлены результаты эффективности распараллеливания алгоритма.

В §1.1 рассматривается постановка задачи для класса решаемых электромагнитных задач. При разработке комплекса программ ЕМ\¥8о1уегЗБ ис-

8

ходной постановкой задачи было найти решение системы уравнений Максвелла в оптической задаче в отсутствии внешних зарядов:

Г _ <ЭВ

гоЖ = —— - Лм - он В от

гоШ. = ^ + J - (ТдБ от

(ИуЪ = 0 ,, ч

V1;

¿гьО = 0 В = рН

>,

где е = е(г, и), ¡1 = ¿¿(г, и) — диэлектрическая и магнитная проницаемости, зависящие от частоты с электрическими потерями, а .7 и Зм — плотность электрического тока и плотность тока магнитного заряда соответственно. Так же будем считать что, материалы могут иметь магнитные и электрические потери энергии, определяемые не зависящими от частоты значениями электрической и магнитной проводимостей &е и <уц■

Комплекс позволяет рассматривать постановки задач, где граничные условия и параметры среды могут иметь вид представленный в Таблице 1. Реализуемая функциональность комплекса программ позволяет моделировать широкий круг задач, таких как рассеивание света на нанообъектах в открытой области, распространение электромагнитных волн в волноводе с диэлектрическими вставками, рассеивание света на фотонных кристаллах и задач, связанных с материалами, имеющими отрицательный коэффициент преломления. Например, задачу моделирование рассеивания волн в открытых микрорезонаторах на основе слоистых метаматериалов [1], а также задача прохождения электромагнитной волны через идеальную линзу.

Граничны е условия

Начальные условия

Модели сред

Условия первого рода.

Задание начального профиля электрического и магнитного поля.

Любая компонента электромагнитного паля равняется заданной функции l(i,y.z.()-J<t,y.lrt г ев #(i,j\i,l)-g(*,J>.i,0; tàtG-zpaxuqao&iacmii

Ê(jr,j',z,f = 0) = /(*,.>>,z); H(x,y,z,t = 0) = g(x,y,z);

Периодические условия.

Метод пвпнвго/расевяннвго поля.

Решение задачи может искаться в виде суммы

Епояное -Епадзющее+Ерассеяниое,

что позволяет моделировать

падение из бесконечности

плоскопараллельной

волны под произвольным заданным

углом._

Г = £(?,«),/!

UPML поглощающий слой, позволяющий решать задачи в неограниченной области.

Поглощающий слои

Дисперсная модель Друде. Позволяет моделировать метаматериалы.

Таблица 1. Функциональность комплекса программ ЕМ\\!Бо 1уе.гЗИ

В §1.2 рассматриваются модели диэлектрической среды для моделирования метаматериалов. Для описания дисперсных сред и метаматериалов предложено использовать модель Друде для е(ш) и ¡л(и>), имеющую вид [9]:

где и ре, Шрт — электрическая и магнитная плазменные частоты, а 7е,7т — частоты столкновений частиц, характеризующие электромагнитные потери в материале.

В §1.3 рассматриваются различные способы описания граничных условий для нестационарных задач электродинамики в ограниченной и неограниченной областях. Для задач в неограниченной области рассматриваются различные формулировки поглощающего слоя РМЬ.

В §1.4 рассматриваются возможные методы задания источника электромагнитных волн при численном решении полной системы уравнений Максвелла с помощью РБТО-метода.

В §1.5 описан метод конечных разностей во временной области численного решения полной системы уравнений Максвелла. Приводится описание расположения электромагнитных полей и контуров в ячейке Йи. Используемая аппроксимация не требует специальных условий сшивки на границах раздела, поскольку выполнение этих условий следует из уравнений Максвелла в интегральной форме. Разностная дивергенция роторов полей Е и Я сохраняется точно. Для описанного метода получена явная РБТБ-схема, которая была использована при создании комплекса программ численного решения трёхмерной системы уравнений Максвелла ЕМ\У8о1уегЗВ.

В §1.6 исследуются основные свойства и параметры 11РМЬ слоя, влияющие на эффективность поглощения. С помощью комплекса программ ЕМ\УБо1уегЗВ проведен анализ параметров 11РМЬ слоя. Показано, что для класса решаемых задач, для обеспечения требуемой точности, вполне достаточно иРМЬ слоя размера р = 20 ячеек при параметрах функции профиля затухания о™\х) = {{х - = -(т+ЭДД(0)1 при т = 4,

1п{В.(0)) = 16 и отдалённости слоя от объектов рассеивания на расстояние порядка длины волны.

В §1.7 приводится описание метода полного/рассеянного поля, используемого в разрабатываемом комплексе программ для задания падающей волны. Данный метод состоит в разделении расчётной области на две подобласти: область полного поля и область рассеянного поля. На границах подобластей ставятся специальные условия сшивки.

В §1.8 приведена численная реализация модели Друде для метаматериа-лов в РОТБ-схеме с учётом разрывных диэлектрической и магнитной прони-цаемостей. Для реализации дисперсного случая в РВТБ-методе используется схема с введением дополнительного уравнения, которая приводится в работе [9].

В §1.9 предложен подход для разработки параллельного алгоритма чис-

ленного решения уравнений Максвелла FDTD-методом, основанный на декомпозиции данных. Параллельный алгоритм для FDTD-метода разработан для суперкомпьютеров с топологией коммуникационной сети 3D-Top или 3D-mesh (в частности IBM BlueGene /Р) и основан на параллелизме по данным. Подход заключается в том, что одна операция выполняется сразу над всеми элементами массива данных. Различные фрагменты такого массива обрабатываются на разных узлах параллельной вычислительной системы. FDTD-метод имеет высокий потенциал к распараллеливанию при декомпозиции расчётного объёма между вычислительными узлами системы, обусловленный сильной локальностью производимых вычислений [20], и позволяет оперировать тера-байтными данными. Исследованы различные способы декомпозиции данных, и проанализирована их эффективность. Рассмотрены возможные варианты балансировки между вычислениями и пересылками данных между узлами вычислительной системы.

В §1.10 рассматриваются технологии программирования при гибридном подходе к программной реализации параллельного алгоритма. В качестве технологии параллельного программирования был применен гибридный подход MPI/OpenMP.

В §1.11 предложен параллельный FDTD алгоритм при 1D декомпозиции для трехмерных задач с протяженной областью. Преимуществом такой схемы является достаточная эффективность и масштабируемость алгоритма без привязки к конкретной топологии коммуникационной сети массивно-параллельного вычислительного комплекса. Однако, такое разбиение возможно применить только к областям протяженным вдоль оси z, иначе на один узел вычислительной системы будет приходиться подобласть, размер которой по z менее 3-х двухмерных слоев, и задача сведется к пересылке граничных значений.

В §1.12 предложен параллельный FDTD алгоритм при 3D декомпози-

ции для произвольных трехмерных задач. Главным преимуществом такой схемы является лучшая по сравнению с 1D декомпозицией эффективность и масштабируемость алгоритма на массивно-параллельных вычислительных системах, имеющих топологию коммуникационной сети 3D-Top или 3D-Mesh. Тем не менее, увеличение количества пересылок и уменьшение объема пересылаемых данных влечет за собой усложнение задачи балансировки между вычислениями и асинхронными пересылками. Если производимый в параллельных ОрепМР циклах расчет занимает время большее, чем пересылка и прием данных с граничных подобластей, то параллельный FDTD алгоритм выполняется эффективно.

В §1.13 рассматриваются особенности целевой архитектуры вычислительной системы IBM Blue Gene /Р, влияющие на эффективность параллельного FDTD алгоритма.

В §1.14 показаны результаты эффективности распараллеливания разработанных алгоритмов численного решения. Разработанные параллельные алгоритмы использовались в комплексе программ EMWSolver3D. Исходный код, а также инструкция по использованию расположены в открытом репо-зитории

по адресу https://github.com/AlexeySemenov/EMWSoIver3D.

В §1.14.1 исследуется эффективность распараллеливания гибридной ОрепМР части. Полученный график ускорения — линейный, и при увеличении числа нитей ускорение максимально. В §1.14.2 исследуется ускорение алгоритма. Характеристикой эффективности распараллеливания кода при решении задачи фиксированного размера является ускорение, которое отражает масштабируемость в так называемом сильном смысле S и может быть определено как [23]:

т(ыприц.,му к '

где , М) — время решения задачи фиксированного размера с вычис-

лительными затратами М на некотором числе процессоров , которое

принято за базовое, а Т(Л>фот(., М) — время решения той же задачи на Атпрщ, процессорах. Исследовано, во сколько раз можно ускорить решение задачи с помощью РБТБ алгоритма при увеличении числа используемых процессорных узлов от 1 до 512 для Ш декомпозиции и от 32 до 2048 для ЗБ декомпозиции. Для Ш декомпозиции расчет проводился для задачи с сеткой 128 х 128 х 1024. С таким объемом данных задача может полностью поместиться на одном вычислительном узле с 2СЬ оперативной памяти. Начиная с 64 узлов ускорение начинает уменьшаться по сравнению с линейным. Однако, при переходе с 256 на 512 узлов ускорение относительно 256 узлов снова линейно. Падение ускорения вызвано тем, что из-за небольшого исходного размера задачи, при разбиении задачи по узлам, нагрузка на каждый вычислительный узел становится крайне мала. Для задачи с сеткой 1024 х 1024 х 4096 ускорение при переходе с 64 на 128 и с 128 на 256 - линейно. Исследована сильная масштабируемость для ЗБ декомпозиции. Ускорение при двукратном увеличении числа узлов варьируется в пределах 1.6-1.8 раз (см. Рис. 1(а)). В §1.14.3 исследована слабая масштабируемость алгоритма. Найден максимальный размер задач, которые можно решать на доступных процессорных мощностях, и получена оценка эффективности использования вычислительных ресурсов при увеличении размера задачи. На Рис. 1(6) представлены графики меры слабой масштабируемости для серии расчётов задач с разностной сеткой; 1024 х 1024 х 2048 и 1024 х 1024 х 8196.

N проц. узлов

1-1 1 и • 1

$1 к ¡за*** -

& 1.03 , Э 0,95 • N Г- . г®»

1 8 2 £ 5 2 10 24 20 «8

N проц. узлов

Рис. 1(а). Ускорение РВТВ-алгоритма при различном числе узлов на сетке 1024 х 1024 х 4096 в случае ЗВ декомпозиции Рис. 1(6) Масштабируемость РВТВ-алгоритма на различных размерах сеток в случае ЗВ декомпозиции

Как видно из графика на Рис. 1(6), параллельный алгоритм при использовании ЗБ декомпозиции хорошо масштабируется и может быть использован для решения электродинамических задач на очень больших сетках.

Результаты первой главы опубликованы в работах [2, 4, 7, 9].

Во второй главе диссертации рассматривается применение комплекса программ для решения задач электродинамики с использованием метамате-риалов.

В §2.1 рассматриваются примеры численного решения задач, позволяющие верифицировать правильность работы компонент программного комплекса. В §2.1.1 рассматривается задача с точечным источником. Точное решение задачи о точечном источнике для Н, компоненты представляет из себя функцию Ханкеля и имеет следующий вид:

Нф,у) =

где 1т - постоянный ток, /3 = ы^/Це, р - расстояние до источника. Решение двухмерной численной задачи для ТЕ-моды рассчитывалось на сетке ЮЛ

х ЮЛ при Ах = Ду = А/10. Источник с частотой бвНг располагался в центре расчетной области х = 5Л, у — 5Л. Было вычислено 200г периодов по времени, чтобы добиться установившегося состояния. При рассмотрении нормированного решения погрешность вычислений на всей рассматриваемой области не превышала 0,00004, что соответствует 0{Ь? +т2). Как следствие поглощения рассеянной волны слоем иРМЬ, видно отсутствие отраженных от границы волн. В §2.1.2 проводилось моделирование прохождения плоской волны через плоскопараллельный слой мстаматериала. Проведено сравнение с аналитическими результатами. Значения коэффициента преломления сред, полученные численного решения, отличаются от теоретического значения не более чем на 0.000032.

Параграф §2.2 посвящён исследованию свойств открытого микрорезонатора Фабри-Перо, содержащего слой метаматериала с отрицательным показателем преломления. В §2.2.1 представлена численная модель микрорезонатора. На основе разработанного комплекса программ проведено численное моделирование открытого резонатора с метаматериалом при помощи разностной РБТБ схемы с учетом частотной дисперсии и поглощающих граничных условий иРМЬ. Аналогичный резонатор рассматривался ранее, однако его размеры существенно превышали длину волны, поэтому представлялось интересным проверить, обладают ли микрорезонаторы теми же уникальными свойствами, что и описанные в работах [21], [22]. В §2.2.2 описываются электромагнитные свойства резонатора с метаматериалом. Численный эксперимент подтверждает предположение о том, что использование метаматериалов позволяет увеличить добротность микрорезонатора на несколько порядков.

В §2.3 рассматривается пример численного моделирования идеальной линзы [9] с помощью программного комплекса ЕМ\У8о1\'егЗБ. Исследована зависимость субволновой разрешающей способности идеальной линзы из метаматериала от её физических параметров. Найдены значения параметров,

при которых возможно получить субволновую разрешающую способность.

Результаты второй главы опубликованы в работах [1-3].

В третьей главе диссертации рассматривается применение программного комплекса для решения задачи нахождения коэффициента прохождения внутри прямоугольного волновода.

В §3.1 рассматривается постановка задачи рассеивания электромагнитной волны на диэлектрическом объекте в волноводе. Требуется найти значение коэффициента прохождения Р электромагнитных волн для различных положений неоднородности, значений диэлектрической проницаемости вставки и значений неоднородности е(х). Коэффициент F находится как отношение комплекснозначных амплитуд волновой моды НЮ до прохождения вставки и после. Для этого рассмотрена задача нахождения амплитуды ведущей волновой моды НЮ из полноволнового РБТБ решения, полученного в ходе численного эксперимента с помощью комплекса программ ЕМ\¥8о1уегЗО.

В §3.2 рассматриваются методы нахождения коэффициента прохождения Р из численного полноволнового решения. Получены уравнения для аналитического решения разностной задачи в случае распространения волновой моды, необходимые для нахождения амплитуды НЮ волновой моды.

В §3.3 произведён расчёт коэффициента прохождения для однородной вставки. Рассмотрена задача для вставки толщины 0.2А и е = 1.9. В соответствии с уравнением, представленным в §3.2 для коэффициента Р в случае однородной вставки, аналитическое значение Не Р = 0.904, а значение, полученное с помощью предложенного метода составляет ЯеР = 0.902. Для той же задачи при е = 2.9 аналитическое значение коэффициента распространения Де F = 0.631, в то время как полученное численно значение Де F = 0.675.

В §3.4 исследована зависимость коэффициента прохождения от физических параметров неоднородности внутри вставки. Рассматривается задача, где диэлектрическая проницаемость е в областях С} п С>1 различна (см. Рис.

2). Толщина диэлектрического слоя <5 равна Л.

Рис. 2: Геометрия и расположение неоднородности (51 внутри диэлектрической вставки <5.

Рис. 3: График коэффициента Г на комплексной плоскости для трёх положений в, = 22 — го = А/12, Л/2, А неоднородности <31 относительно оси г при £2 = 1.9..15.9 с шагом 1

Размеры неоднородности ж2 - х\ = уг — = ¿2 — г\ — А/10. Расчёты проводились для трёх положений й = — ¿о = А/12, А/2, А неоднородности <51 относительно оси г. Значение диэлектрической проницаемости для вставки <5 было выбрано следующим: £1 = 1,9, а значение диэлектрической проницаемости £2 для неоднородности менялось в каждом расчете от 1.9 до 15.9 с шагом 1. Для каждого расчёта был посчитан коэффициент Р. На Рис. 3 представлен график коэффициента ^ на комплексной плоскости. Каждая точка графиков на комплексной плоскости - полученное численное значение ^ для фиксированных параметров неоднородности <1 и £2- Кривая образуется из точек для фиксированного положения неоднородности <1 и последовательно увеличивающегося значения параметра е^, меняющегося в пределах от 1.9 до 15.9 с шагом 1. Отсутствие самопересечений кривых на графике говорит о возможности однозначного определения диэлектрической проницаемости

£2 по коэффициенту прохождения F. Полученные значения коэффициента F можно использовать при решении обратной задачи определения диэлектрической проницаемости, если необходимо многократное решение прямой задачи нахождения F.

Результаты третьей главы опубликованы в работах [5, 6, 8].

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

Основные результаты

• Разработан программный алгоритм для решения электродинамических задач, использующий гибридные технологии MPI и ОрепМР. Создан многофункциональный комплекс программ на многоядерных вычислительных системах и суперкомпьютере IBM BlueGene/P, позволяющий решать широкий круг электродинамических задач. Показано, что созданный комплекс имеет хорошую масштабируемость и может эффективно выполняться на массивно-параллельных системах, оперируя с терабайтами вычислительных данных.

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

• С помощью комплекса впервые решена задача вычисления добротности открытого резонатора Фабри-Перо со вставкой из метаматериала, ранее проводившегося только для не субволновых размеров резонатора.

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

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

Список публикаций

1. Семенов А. Н., Смирнов А.П., Игнатьева Д. О., Сухоруков А.П. Математическое моделирование открытого микрорезонатора со слоем метамате-риала // Известия РАН. Серия физическая 75(12), 2011, с. 163 7-1640

2. Семенов А.Н., Смирнов А.П. Численное моделирование уравнений Максвелла с дисперсными материалами //Математическое моделирование, 2013, том 25, номер 12, с. 19-32

3. Smirnov А.Р., Semenov A.N., Ignatyeva, D.O., Sukhorukov. A.P. Full-wave modeling of open subw avelength resonator with metamaterial // Progress in Electromagnetics Research Symposium , Moscow, 2012, pp. 1254-1258

4. Smirnov, A.P., Semenov, A.N. Full wave Maxwell's equations solver EMWSolver3D //' Progress in Electromagnetics Research Symposium, Moscow, 2012, pp. 252-255

5. Smirnov A. P., Semenov A.N., Shestopalov Y. Modeling of Electromagnetic Wave Propagation in Guides with Inhomogeneous Dielectric Inclusions // Inverse problems and large-scale computations/ Eds. Larisa Beilina, Yury V. Shestopalov, Springer

6. Smirnov A.P., Semenov A.N., Shestopalov Y. FDTD Simulation of Waveguide

with Non-uniform Dielectric Slab // Progress in Electromagnetics Research Symposium , Stockholm, 2013, pp. 76-78

7. Семенов A. H. Параллельная реализация численного решения уравнений Максвелла Fdtd методом для больших задач // Научный сервис в сети Интернет: экзафлопсное будущее: Труды Международной суперкомпыо-терной конференции М.: Изд-во МГУ, 2011. с.530-533

8. Smirnov А.P., Semenov A.N., Shestopalov Y. Modeling of electromagnetic wave propagation in guides with inhomogeneous dielectric inclusions // The Workshop proceedings (tentative title "Inverse Problems and Large-Scale Modeling with Applications in Electromagnetics. Proceedings of conferences supported by the Visby Program"), Karlstad, Sweden, 2012, pp. 13-14

9. Семенов A.H., Смирнов А.П. Комплекс программ численного решения уравнений Максвелла EMWSOLVER3D //VI Всероссийская конференция "Радиолокация и радиосвязь" , JRE (Журнал Радиоэлектроники) -ИРЭ им. В.А.Котельникова РАН Москва, 19-22 ноября 2012 г, с. 146-150

Цитированная литература

1. A.Tafiove and S.C.Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method // Norwood, MA:Artech, 2000.

2. Z.S.Sacks, D.M.Kingsland, R.Lee, and J.F.Lee, A perfectly matched anisotropic absorber for use as an absorbing boundary condition // IEEE Trans.Antennas and Propagation, 1995, vol.43 ,no,12, pp.1460-1463.

3. Kane Yee, Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // Antennas and Propagation, 1966, IEEE Transactions on 14, pp. 302-307

4. S.D. Gedney, An anisotropic perfectly matched layer absorbing media for the truncation of FDTD latices // Antennas and Propagation, 1996, IEEE Transactions on 44, pp.1630-1639.

5. A. Taflove, J. Dabkowski, and M. Genge, Mitigation of buried pipeline voltages due to 60 Hz ac inductive coupling. Part I: Design of joint rights of way; Part II: Pipeline grounding methods // IEEE Trans. Power Apparatus Systems, 1979, vol. 98, pp. 1806-1823.

6. А. Зеленин, Решение уравнений Максвелла методом FDTD // По материалам сайта http://zfdtd.narod.ru

7. А.Н. Боголюбов, И.А. Буткарев, Ю.С. Дементьева Численное моделирование двумерных фотонных кристаллов // "ЖУРНАЛ РАДИОЭЛЕКТРОНИКИ'^ 11,2006

8. David М. Hockanson Perfectly Matched Layers Used as Absorbing Boundaries in a Three-dimensional FDTD Code

9. Yang Hao, Raj Mittra FDTD Modeling of Metamaterials Theory and Applications // Artech House Publishers, 2009.

10. Simpson. 11. and A Taflove Three-dimensional FDTD modeling of impulsive ELF propagation about the Eanh-sphere // IEEE Trans. Antennas Propagai, Vol. 52, 2004. pp. 443-451

11. Simpson. J J. and A. Taflove Efficient modeling of impulsive ELF antipodal propagation about the Earth-sphere using an optimized two-dimensional geodesic FDTD grid // IEEE Antennas Wireless Propagai Leu.. Vol. 3. 2004. pp. 215-218.

12. Hagness. S C. A. Taflove. and J. E Bridges. Three-dimensional FDTD analysis of a pulsed microwave cocfocal system for breast cancer detection design of an antenna-array element // IEEE Trans Antennas Propagaition Vol 47. 1999. pp. 783-791.

13. Bond. E J. X Li. S. C. Hagness. and B D Van Veen Microwave imaging via ipace time beamforming for early detection of breast cancer // IEEE Trans Antennas Propagat. Vol 51. 2003. pp 1690-1705

14. Simpson. J. J.. A Taflove. i A Mu. and H. Heck Computational and experimental study of a microwave electromagnetic bandgap structure with waveguiding defect for potential use as a bandpass wireless interconnect // IEEE Microwave Wireless Components Lett. Vol 14. 2001. pp 345-345

15. Simpson. J. J.. A Taflove. J. A. Mix. and H Heck Advances in hypcrspecd digital interconnects using electromagnetic bandgap technology: Measured low-loss 43-GHz passband centered at 50 GHz // Proc IEEE Antennas Propagai. Soc Inti Svmp . Washington. D.C, 2005

16. Park. H.C., S.H. Kim. S.H Kwon. Y.O. Ju. J.K. Yang. J.H. Back. S.B Kim. and Y -II Lee nElectrically driven tingle-cell photonic crystal laser // Science, Vd 305. 2004. pp. 1444-1447

17. Yanik. M F . S Fan. M. Soljactc. and J. D. Joannopoulos All-optical transistor action with bistable switching in a photonic crystal cross-waveguide geometry // Optics Lett., Vol. 28. 2003. pp 2506-2508.

18. Ardavan F. Oskooi, David Roundy, Mihai Ibanescu, Peter Bermel, J. D. Joannopoulos, and Steven G. Johnson MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method // Computer Physics Communications 181, 2010, pp.687-702 .

19. А. V. Zakirov, V. D. Levchenko CFmaxwell - The Program Code for FDTD Modelling of Very Large Size Problems //in Progress in Electromagnetics Research Symposium, Stockholm, Sweden, August 2013

20. Wenhua Yu, Raj Mittra, Tao Su, Yongjun Liu, Xiaoling Yang Parallel Finite-Difference Time-Domain Method // Artech House; 1 edition, 2006

21. Д. О. Сапарина, А. П. Сухоруков // Известия вузов. Прикладная нелинейная динамика. 2009, Т. 17, поЗ. С. 3-16

22. Saparina, D. О.; Sukhorukov, А. Р. // Laser Physics. 2009. Vol. 19. No. 5. pp.1125-1130

23. P. Chatelain et al. Billion vortex particle direct numerical simulations of aircraft wakes» // Comput. Methods Appl. Mech. Engrg., 2008, Vol.197, pp. 1296-1304.

Напечатано с готового оригинал-макета

Подписано в печать 21.11.2013 г. Формат 60x90 1/16. Усл.печ.л. 1,0. Тираж 80 экз. Заказ 398.

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к. Тел. 8(495)939-3890/91. Тел./факс 8(495)939-3891.

I 1

Московский государственный университет им. М.В. Ломоносова. Факультет вычислительной математики и кибернетики

На правах рукописи

04201452320

Семенов Алексей Николаевич

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

систем

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель

к. ф.-м. II., доц.

Смирнов Александр Павлович

Москва - 2013

Содержание

Введение 4

Глава 1. Создание комплекса программ численного моделирования уравнений Максвелла 10

1.1. Класс решаемых электромагнитных задач ............10

1.2. Модели диэлектрической среды...................11

1.3. Граиичные условия..........................12

1.4. Источники электромагнитных волн ................14

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

1.6. Использование поглощающего слоя при моделировании открытых задач...............................24

1.7. Метод нолного/рассеяного поля ..................30

1.8. Численная реализация модели Друде для метаматериалов в FDTD схеме..............................31

1.9. Параллельный алгоритм численного решения уравнений Максвелла FDTD-методом ........................33

1.10. Гибридный подход к программной реализации параллельного алгоритма...............................34

1.11. 1D декомпозиция для трехмерных задач с протяженной областью 35

1.12. 3D декомпозиция для произвольных трехмерных задач.....38

1.13. Особенности целевой архитектуры вычислительной системы . . 41

1.14. Эффективность распараллеливания разработанных алгоритмов численного решения.......................45

Глава 2. Применение комплекса программ для решения оптических задач с использованием метаматериалов 54

2.1. Апробация в базовых оптических задачах.............55

2.2. Моделирование открытого микрорезонатора со слоем метама-териала................................62

2.3. Моделирование идеальной линзы..................73

Глава 3. Применение комплекса программ для решения задачи оценки коэффициента прохождения волноводной моды

в прямоугольном волноводе 78

3.1. Моделирование электромагнитного волновода с диэлектрической вставкой.............................78

3.2. Получение коэффициента прохождения из полноволнового решения .................................81

3.3. Расчёт коэффициента прохождения для однородной тонкой вставки ................................86

3.4. Исследование зависимости коэффициента прохождения от физических параметров неоднородности внутри вставки......90

Заключение 92

Литература 93

Введение

Актуальность работы. При моделировании распространения электромагнитных волн в радиотехнике [22],[23],[8], микроэлектронике [28],[27], задачах нанооптики [29],[30], [10] и биодиагностики [25], [26], требуется решение полной трёхмерной системы уравнений Максвелла на сетках больших размеров. При реализации на вычислительной системе размер таких сеток составляет порядка терабайта и более. Большие сетки необходимы для описания исоднородностей, характерный размер которых значительно меньше длины волны, а общий размер исследуемой системы равен десяткам длин волн. Для решения таких задач необходимы вычислительные системы с большим объёмом памяти и высокой производительностью: суперкомпьютеры с массивно-параллельной архитектурой. Для их использования требуются специализированные комплексы программ, способные эффективно использовать терабайтные объёмы оперативной памяти и более чем тысячи процессоров вычислительной системы. При реализации комплекса программ необходимо также учитывать специфику архитектуры конкретной вычислительной системы. Комплекс должен быть многофункциональным и рассчитанным на решение задач как в ограниченной, так и в не ограниченной области.

Определение электромагнитных полей в произвольной, неоднородной диэлектрической среде - это важный практический предмет исследования волновых эффектов, в том числе и на микро и нано уровне. Аналитические решения получены, как правило, для объектов простой формы. Для объектов сложной структуры необходимо прибегнуть к численному моделированию задачи. Численное решение полной системы уравнений Максвелла позволяет получить значения компонент электромагнитного поля для последующего исследования и анализа электромагнитных свойств моделируемых структур.

Одним из наиболее распространённых способов численного решения пол-

ной системы уравнений Максвелла является явный метод конечных разностей во временной области, в иностранной литературе fmitc-difFcrencc timedomain (FDTD) [1]. Метод основан на использовании сетки Йи [4], в узлах которой располагаются компоненты электромагнитного поля. Основываясь на уравнениях Максвелла в интегральной форме, можно построить их конечно-разностную аппроксимацию. Поскольку FDTD-мстод решает задачу во временной области, он позволяет решать импульсные и негармонические задачи и получить результат для широкого спектра длин волн за один расчёт. Это необходимо при решении задач, в которых неизвестны резонансные частоты, или в случае моделирования широкополосных сигналов. Также FDTD-метод позволяет проследить временную эволюцию распространения волны в моделируемом объёме, что может быть необходимо для детального исследования процесса формирования устойчивой картины распределения электромагнитных волн. Метод позволяет непосредственно моделировать краевые эффекты и эффекты экранирования и удобен при задании анизотропных, дисперсных и нелинейных сред.

При численном решении FDTD-методом задач в неограниченном пространстве на границах расчётной области чаще всего используется искусственный поглощающий слой PML (Perfectly matched layer) [11], [9]. Рассеянная волна, уходя на бесконечность, поглощается анизотропным слоем на границе и, при правильном выборе параметров поглощающего слоя и разностной сетки, отражение от PML слоя минимально [5],[3].

Величина шага дискретизации по пространству должна быть значительно меньше исследуемых длин волн и размера исследуемой структуры. Это может потребовать сеток с маленьким шагом по пространству, что означает большие затраты памяти и возросшее время расчёта. При моделировании дисперсных и анизатропных сред также растут необходимые вычислительные ресурсы. Актуальной становится задача создания программного комплекса

на высокопроизводительных массивно-параллельных вычислительных системах с большим объёмом оперативной памяти, который позволит эффективно решать класс электродинамичских задач, описываемых полной системой уравнений Максвелла.

Существует ряд комплексов, реализующих параллельный FDTD-мстод на кластерных вычислительных системах([3б], [37]). Однако, использование таких комплексов, как правило, рассчитано на массивно-параллельные вычислительные системы общей архитектуры, что делает невозможным эффективно использовать их на специализированных суперкомпьютерах и полностью задействовать их вычислительные мощности. На факультете ВМК Московского государственного университета имени М.В. Ломоносова установлен суперкомпьютер IBM BlueGene/P, имеющий 8192 процессоров и 4 терабайта оперативной памяти. Разработанный в диссертации программный комплекс был эффективно реализован для этого суперкомпьютера и вычислительных систем схожей архитектуры.

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

Научная новизна. Разработан программный алгоритм для решения крупномасштабных электродинамических задач, который с помощью гибридных технологий MPI и ОрепМР может эффективно выполнятся на современных массивно-параллельиых вычислительных системах. На основе предложенного алгоритма создан программный комплекс на многоядерных вычислительных системах и суперкомпьютерах серии IBM Blue Gene/P. Показано, что комплекс имеет хорошую масштабируемость и может эффективно

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

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

На защиту выносятся следующие основные результаты и положения:

• Разработан программный алгоритм для решения электродинамических задач, использующий гибридные технологии MPI и ОрепМР. Создан многофункциональный комплекс программ на многоядерных вычислительных системах и суперкомпьютере IBM BlueGene/P, позволяющий решать широкий круг электродинамических задач. Показано, что созданный комплекс имеет хорошую масштабируемость и может эффективно выполняться на массивно-параллельных системах, оперируя с терабайтами вычислительных данных.

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

• С помощью комплекса впервые решена задача вычисления добротности

открытого резонатора Фабри-Перо со вставкой из метаматериала, ранее проводившегося только для не субволновых размеров резонатора.

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

Апробация работы. Основные результаты диссертации докладывались на:

• Международной конференции "Progress in Electromagnetics Research Symposium"(Stockholm, Sweden, 2013 год)

• Международной конференции "Progress in Electromagnetics Research Symposium" (Moscow, Russia, 2012 год)

• Международном научном семинаре "Inverse Problems and Large-Scale Modeling with Applications in Electromagnetics" (Karlstad, Sweden, 2012 год)

• VI Всероссийской конференция "Радиолокация и радиосвязь" Российская академия наук Институт радиотехники и электроники им. В. А. Котелышкова (Москва, 2012 год)

• Объединенном Фельдовском семинаре по электродинамике и антеннам. (Москва, 2012 год) Российская академия наук Институт радиотехники и электроники им. В. А. Котелышкова

• XIII Всероссийской школы-семинара "Физика и применение микроволн" (Звенигород, МГУ им. М.В. Ломоносова, 2011 год)

• Международной суперкомпьютерной конференции "Научный сервис в сети Интернет: экзафлопеное будущее" (Новороссийск, 2011 год)

• Научной конференции "Ломоносовские чтения"(Москва, МГУ им. М.В. Ломоносова, 2011 год)

Публикации. Материалы диссертации опубликованы в 9 печатных работах, из них б статей в журналах списка ВАК ([40]-[45]) и 3 тезиса докладов ([46]-[48]).

Структура и объем диссертации Диссертация состоит из введения, 3 глав, заключения и библиографии. Общий объем диссертации 100 страниц, из них 89 страниц текста, включая 23 рисунка. Библиография включает 60 наименований на 7 страницах.

Глава 1

Создание комплекса программ численного моделирования уравнений Максвелла

В диссертационной работе автором был создан комплекс программ ЕМЛУБокегЗБ, предназначенный для численного решения полной системы уравнений Максвелла.

1.1. Класс решаемых электромагнитных задач

При разработке комплекса программ ЕМ\¥8о1уегЗО исходной постановкой задачи было найти решение системы уравнений Максвелла в оптической задаче в отсутствии внешних зарядов:

гоШ

дВ

Зм — 0яВ

ал)

дг

гоШ = —- + 3 - (7ЕВ дЬ

(НуВ = О (НуВ = О В = цН Б = еЕ

где £ = е(г, а;), /2 = /г(г,а;) — диэлектрическая и магнитная проницаемости, зависящие от частоты с электрическими потерями, а J и Ям ~ плотность электрического тока и плотность тока магнитного заряда соответственно. Так же будем считать что, материалы могут иметь магнитные и электрические потери энергии, определяемые не зависящими от частоты значениями электрической и магнитной проводимостей ое и Ое-

Граничные условия Начальные условия Модели сред

Условия первого рода. Любая компонента электромагнитного поля равняется заданной функции Ё{х,у,г,1)-Дх,у,2Л геС Й(х,у,2,1)-гдеС-гратцаобластц Задание начального профиля электрического и магнитного ПОЛЯ. Е(х,у, /,Г = 0) = /(х,>%2); Н(х;у,гЛ = 0) = §(х:у}г); е=е(г), М=Р(Г)

£ = = Ц(г,ь») Дисперсная модель Друде. Позволяет моделировать метаматериалы.

Периодические условия. Метод полного/рассеянного поля. Решение задачи может искаться в виде суммы Е полное =Епадэющее+Ерассеянное, что позволяет моделировать падение из бесконечности плоскопараллельной волны под произвольным заданным углом.

Поглощающий слой иРМ1. поглощающий слой, позволяющий решать задачи е неограниченной области.

Таблица 1. Функциональность комплекса программ ЕМ\Уво1иегЗВ

Комплекс позволяет рассматривать постановки задач, где граничные условия и параметры среды могут иметь вид, представленный в Таблице 1. Реализуемая функциональность комплекса программ позволяет моделировать широкий круг задач, таких как рассеивание света на нанообъектах в открытой области, распространение электромагнитных волн в волноводе с диэлектрическими вставками, рассеивание света па фотонных кристаллах и задач, связанных с материалами, имеющими отрицательный коэффициент преломления. Например, задачу моделирования рассеивания волн в открытых микрорезонаторах на основе слоистых метаматериалов [40], а также задачу прохождения электромагнитной волны через идеальную линзу.

1.2. Модели диэлектрической среды

При описании свойств материалов с отрицательным коэффициентом преломления нужно иметь в виду, что эти материалы должны обладать частотной дисперсией. Действительно, если е и /х оба отрицательны, то при отсутствии дисперсии плотность энергии электрического и магнитного поля,

равная [15]:

IV =1-{еЕ2 + /х#2), (1.2)

будет отрицательной. При наличии частотной дисперсии плотность энергии (1.2) имеет вид [15]:

дЩи1\Е2 + дЫимн2

2 дш ди

В ряде работ [13] [18] предложено описывать метаматериалы дисперсной моделью Друде [13]. Для модели Друде е{и>) и ¡1{и}) принимают вид:

е(и) = £0 I 1 - ре

О)2 — гШ'Уе

и? \ {1А)

= /¿о ( 1--—- '

•'рт

,2 ^ГТГ

- гш'Ут у

где МрсШрт - электрическая и магнитная плазменные частоты, а 7е,7т - частоты столкновений частиц, характеризующие электромагнитные потери в материале. Заметим, что эта модель удовлетворяет условиям

аКыМ>0а|р(цМ>0- (15)

дсо ди

1.3. Граничные условия

При решении задач с наличием границ используются заданные граничные условия. Например, для описания идеального проводника (РЕС)[1], используется условия равенства нулю тангенциальных компонент электромагнитного поля. Также, в рамках комплекса реализованы периодические граничные условия.

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

используются поглощающие граничные условия. Существует ряд аналитических условий на бесконечности для уравнений Максвелла [31], [32]. Большинство из них имеет узкую область применения и ряд существенных ограничений [34] [33]. При численном решении задач на бесконечности также существуют различные методы вычисления значений поля на границе, зачастую привносящие значительные искажения в полученное решение [1]. Хорошие результаты при численном моделировании различных задач электродинамики показали поглощающие граничные условия Perfect matched layer(PML)[1]. PML требует аккуратного выбора параметров поглощающего слоя и разностной сетки для получения приемлемого результата.

Есть несколько эквивалентных формулировок PML. Первоначальная формулировка J. P. Berenger [2] в зарубежной литературе именуется split-field PML. В ней решение искусственным образом разбивается на сумму двух опять же искусственных компонент электромагнитного поля. На данный момент, более общей формулировкой считается uniaxial PML или UPML [5], в котором в области PML решается тоже волновое уравнение, но в комбинации с искусс