«ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ И КОМПЬЮТЕРНЫЙ ДИЗАЙН ОПТИЧЕСКИХ СВОЙСТВ НАНОСТРУКТУРИРОВАННЫХ МАТЕРИАЛОВ ...»

РОССИЙСКИЙ НАУЧНЫЙ ЦЕНТР

"КУРЧАТОВСКИЙ ИНСТИТУТ"

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

Дейнега Алексей Вадимович

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ И КОМПЬЮТЕРНЫЙ

ДИЗАЙН ОПТИЧЕСКИХ СВОЙСТВ

НАНОСТРУКТУРИРОВАННЫХ МАТЕРИАЛОВ

Специальность 01.04.05 оптика

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

Научный руководитель кандидат физико-математических наук доцент Потапкин Б. В.

Москва 2010 2 Оглавление Введение 1 Развитие метода решения уравнений Максвелла в конечных разностях 1.1 Обзор литературы........................ 1.2 Метод подсеточного сглаживания............... 1.3 Итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру....... 1.4 Выводы.............................. 2 Применение шаблонного метапрограммирования для эффективной реализация FDTD 2.1 Контурный подход к дискретизации уравнений Максвелла. 2.2 Код как совокупность взаимодействующих компонент.... 2.3 Стадии расчета.......................... 2.4 Оптимизация использования памяти.............. 2.5 Проведение параллельных расчетов.............. 3 Металлические фотонные кристаллы как источники светового излучения 3.1 Обзор литературы........................ 2.6 Выводы.............................. 3.2 Оптимизация излучательных характеристик металлических фотонных кристаллов в видимом диапазоне.......... 3.3 Учет внешней матрицы, необходимой для фиксации фотонного кристалла.......................... 3.4 Выводы.............................. 4 Антиотражающие текстурированные покрытия 4.1 Обзор литературы........................ 4.2 Численное моделирование текстурированных покрытий во всем диапазоне размеров текстуры............... 4.3 Кремниевые текстурированные покрытия........... 4.4 Выводы.............................. 5 Заключение 5.1 Основные результаты и выводы работы............ 5.2 Благодарности.......................... Литература Введение Диссертация посвящена численному моделированию оптических свойств наноструктурированных материалов, а именно, в ней исследуются следующие вопросы:

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

• Оптимизация размеров и формы антиотражающих нанотекстурированных покрытий.

Основным используемым численным методом является метод решения уравнений Максвелла в конечных разностях (Finite-Dierence TimeDomain, FDTD). Для расчетов применяется специально написанная параллельная программа, включающая в себя ряд новых оригинальных численных методов, которые также описываются в данной диссертации.

Актуальность работы. Наноструктурированные материалы представляют собой новый тип материалов, характеризующийся малым размером (характерный размер порядка нескольких сотен нанометров и менее) и сложной организацией составляющих его элементов. Наноматериалы обладают уникальными физическими свойствами, что позволяет им находить применение во многих промышленных областях: в вычислительной технике, энергетике, медицине и т. д. [1] К наноматериалам относятся фотонные кристаллы [2], представляющие собой структуры с периодически меняющейся в пространстве диэлектрической проницаемостью. Оптические свойства фотонных кристаллов сильно отличаются от оптических свойств сплошных сред: электромагнитные волны в фотонном кристалле имеют зонный спектр и координатную зависимость, аналогичную блоховским волнам электронов в обычном кристалле. Перераспределение плотности фотонных состояний в фотонном кристалле приводит к изменению спектра их теплового излучения [3].

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

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

Еще одним примером практического использования наноматериалов являются текстурированные антиотражающие покрытия. Интерес к ним особенно возрос в последнее время, а именно, появилась масса новых работ по успешному изготовлению антиотражающих нанотекстурированных покрытий (см., напр., [4]), используемых, в частности, в солнечной энергетике. Параллельно активно ведется теоретическое исследование оптических свойств текстурированных покрытий (см., напр., [5]). Отметим, что в имеющихся работах, как правило, либо отдельно рассматриваются длинноили коротковолновый предельные случаи, либо видимый диапазон для размеров составляющих текстуру элементов порядка нескольких сотен нанометров. В это же время является актуальным единое исследование текстурированных поверхностей во всем диапазоне длин волн. Такое исследование позволило бы установить, каков должен быть размер составляющих текстуру элементов, чтобы для заданного диапазона длин волн отражение было бы минимальным.

В этой работе основным используемым численным методом является метод решения уравнений Максвелла в конечных разностях (англ. - FiniteDierence Time-Domain method, FDTD) [6]. Выбор этого метода вызван присущей ему высокой параллельной эффективностью, что позволяет применять его для расчета больших задач на кластерных вычислительных системах.

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

Исследование оптических свойств наноматериалов включает в себя получение спектров прохождения и отражения от периодических структур при падении на них плоской волны. Их расчет FDTD для случая нормального падения предполагает применение периодических граничных условий.

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

Цель работы состоит в 1) исследовании излучательных характеристик металлических фотонных кристаллов, а также возможности их практического применения в качестве источников света c высокой селективностью излучения в видимом диапазоне; 2) исследовании антиотражающих свойств текстурированных покрытий во всем диапазоне размеров составляющих текстуру рассеивателей, а также установлении оптимального размера для заданного диапазона длин волн; 3) дальнейшем развитии метода FDTD, включающем, в частности, разработку метода подсеточного сглаживания для дисперсных сред и метода расчета наклонного падения плоской волны на периодическую структуру.

Научная новизна работы.

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

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

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

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

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

5. Предложен новый подход программной реализации метода решения уравнений Максвелла в конечных разностях (FDTD), основанный на интегральном представлении уравнений Максвелла. В рамках этого подхода написана параллельная программа на языке программирования С++.

В ходе ее написания реализован ряд оригинальных решений: оптимизация использования памяти с помощью упаковки используемых в разностных уравнениях коэффициентов; увеличение производительности посредством разбиения алгоритма обновления сеток на основной цикл, поддающийся векторизации, и специализированные поправки к нему; увеличение параллельной эффективности путем балансировки доменов. Специально проведенные тесты продемонстрировали линейный характер масштабируемости программы вплоть до тысячи вычислительных ядер. Программа находится в открытом доступе в интернете (http://fdtd.kintechlab.com).

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

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

Положения, выносимые на защиту.

1. Геометрические параметры вольфрамовых фотонных кристаллов типа прямого опала, при которых КПД источника света на основе фотонных кристаллов выше в несколько раз по сравнению с обычными лампами накаливания.

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

3. Асимптотические зависимости величины отражения от геометрических параметров текстурированной поверхности для длинно- и коротковолнового пределов.

4. Ключевым фактором, влияющим на оптимальный размер составляющих текстуру рассеивателей, является характер замощения подложки их основаниями: в случае полного замощения наименьшее значение отражения достигается при макроскопических размерах рассеивателей, а в случае неполного - оптимальный размер рассеивателей имеет порядок длины волны.

5. Параллельная программа FDTD, которая может использоваться для предсказательного моделирования в разнообразных оптических приложениях.

6. Новые численные методы в рамках FDTD: метод подсеточного сглаживания для дисперсных сред и итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру.

Практическая ценность работы.

1. Найдены геометрические параметры вольфрамовых фотонных кристаллов типа прямого опала, при которых КПД источника света на основе фотонного кристалла оказывается выше в несколько раз по сравнению с обычными лампами накаливания.

2. Установлено, что для создания высокоэффективных источников света на основе фотонного кристалла необходим поиск материалов со слабым поглощением при высоких температурах.

3. Найдена оптимальная геометрия текстурированной поверхности, при которой достигается наименьшая величина отражения для заданного диапазона длин волн.

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

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

6. Предложены и реализованы численные методы, позволяющие существенно увеличить быстродействие и точность расчетов FDTD.

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

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

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

Апробация работы. Основные результаты диссертационной работы были доложены на научно-технических конференциях The International Conference on Computational Science and Its Applications“ (Малайзия, Kuala Lumpur, 2007), Conference on Computational Physics“ (Бразилия, Ouro Preto, 2008), Conference on Computational Physics“ (Тайвань, Kaohsiung, 2009), Многомасштабное моделирование процессов и структур в нанотехнологиях“ (Москва, 2009), VII Курчатовская молодежная научная школа“ (Москва, 2009). 52-ая научная конференция МФТИ (Москва, 2009).

Реализованная программа FDTD стала призером на конкурсе Максимальная масштабируемость“, проведенном в рамках РОСНАНОФОРУМА-2009 компаниями Intel и "Роснано".

Публикации. По материалам диссертации опубликовано 5 работ в реферируемых научных изданиях [8–12], 1 в нерефериуемом научном издании [13] и тезисы российских и международных конференций.

1. Развитие метода решения уравнений Максвелла в 1.1. Обзор литературы Моделирование оптических свойств наноматериалов требует непосредственного численного решения уравнений Максвелла. Одним из наиболее популярных предназначенных для этого методов является метод решения уравнений Максвелла в конечных разностях (англ. - Finite-Dierence Time-Domain method, FDTD) [6]. Базовый алгоритм этого метода был предложен Кейном Йи в 1966 г. [14]. Однако имеющуюся в данный момент популярность FDTD приобрел только в 90х гг. прошлого столетия, когда он стал основным для численного моделирования самых разных оптических приложений. В это же время начало экспоненциально расти и число публикаций, посвященных дальнейшему развитию этого метода.

Существует целый ряд причин, обуславливающих популярность FDTD:

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

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

• Естественно присущая методу высокая параллельная эффективность позволяет применять его для расчета больших задач на кластерных вычислительных системах;

• В FDTD существует возможность как для анализа зависимости оптических характеристик исследуемой структуры от длины волны, так и для моделирования нестационарных процессов;

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

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

В следующем параграфе мы более подробно изложим содержание базового алгоритма FDTD.

1.1.1. Алгоритм Йи В отсутствие свободных зарядов уравнения Максвелла в дифференциальной форме имеют вид:

где E и H - напряженности электрического и магнитного полей, D - электрическая, B - магнитная индукция, J - плотность электрического тока, M - ее магнитный аналог, который мы включили для общности некоторых используемых далее выражений.

В целях простоты изложения мы ограничим рассмотрение только изотропными, недисперсными и линейными среды. Тогда вектора E и D, H и B связаны следующими соотношениями:

где, r, µ, µr - абсолютные и относительные диэлектрическая и магнитная проницаемости среды, 0 и µ0 - электрическая и магнитная постоянные.

Плотность электрического и магнитного токов представим в виде:

где Jsource - плотность тока свободных источников, - электропроводность, Msource и - их магнитные аналоги.

Подставив соотношения (1.5), (1.6) в первую пару уравнений Максвелла (1.1),(1.2), получим Согласно алгоритму Йи [6,14], уравнения Максвелла дискретизуются с использованием центрального разностного отношения для приближения пространственной и временной производных. Для этого сеточные узлы, в которых хранятся компоненты электрического и магнитного полей, смещены по отношению друг к другу на половину сеточного шага по каждой из пространственных переменных (Рис. 1.1). В результате те узлы, которые соответствуют компонентам полей E расположены таким образом, что каждая компонента E окружена четырьмя компонентами H, и наоборот.

Рис. 1.1. Положения сеточных узлов для компонент электрического и магнитного полей на сетке Йи.

Согласно алгоритму Йи соответствующие компонентам E и H узлы сдвинуты относительно друг друга по времени на половину временного шага (в качестве примера это продемонстрировано на Рис. 1.2 для одномерной сетки Йи). Для расчета значений поля E на n + 1/2-ом временном шаге используются значения поля H на n-ом. Аналогичным образом значения поля H на n + 1-ом шаге рассчитываются с использованием значений поля E на n+1/2-ом шаге. Эта процедура продолжается до тех пор, пока расчет не будет закончен.

Прежде чем получить конкретные разностные уравнения, используемые в алгоритме Йи, введем следующие обозначения. Будем ставить в соответствие каждому сеточному узлу три целых числа i, j и k, определяющих положение этого узла в пространстве как где x, y и z суть сеточные шаги по соответствующим направлениям.

Будем обозначать произвольную сеточную функцию u как где t - шаг по времени, а n - текущий временной шаг.

Рис. 1.2. Развертка одномерной сетки Йи по времени: узлы, соответствующие электрическому и магнитному полям, сдвинуты друг относительно по пространству и по времени на половину шага.

Как уже упоминалось, в алгоритме Йи для аппроксимации присутствующих в уравнениях Максвелла производных используется центральное разностное отношение. Так, для производной по координате x в точке (ix, jy, kz) в момент времени nt имеем Аналогично для производной по времени t Применяя эти приближения для производных, участвующих в проекции закона Ампера на ось x получаем i,j+1/2,k+1/ i,j+1/2,k+1/2 i,j+1/2,k+1/2 Ex |i,j+1/2,k+1/2.

Стоящие в правой части (1.14) переменные берутся на временном шаге n, включая поле Ex. Поскольку для значения Ex в момент времени n данных на сетке нет, нужно для него пользоваться каким-либо приближением. Хорошим выбором является усреднение по соседним временным слоям:

Положим = x = y = z. Тогда после подстановки (1.15) в (1.14), мы можем явно выразить Ex на n + 1/2 шагу:

Ex |i,j+1/2,k+1/2 = Ca,Ex |i,j+1/2,k+1/2 Ex |i,j+1/2,k+1/ где Мы получили разностное уравнение, соответствующее проекции закона Ампера на ось x, которое вместе с пятью оставшимися аналогичными разностными уравнениями и формируют алгоритм Йи.

Выбор значений x, y, z обуславливается геометрией задачи и спектральным составом излучения. Можно дать следующую рекомендацию: на характерный размер объекта (радиус шариков, толщину экрана и т. п.) должно приходится не менее нескольких сеточных шагов, а на характерную длину волны - от десяти и больше. Величина значения t ограничена сверху условием Куранта выполнение которого необходимо для устойчивости разностной схемы.

В этом параграфе мы рассмотрели реализацию алгоритма Йи для случая, когда среда характеризуется не зависящими от частоты скалярными значениями, µ,,. Численное моделирование анизотропных, дисперсных и нелинейных сред требует модификации этого алгоритма и применения вспомогательных разностных уравнений [6, 15, 16].

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

свойств наноматериалов Перечислим основных участников“ численного эксперимента FDTD.

Материальные тела, оптические свойства которых мы исследуем.

Если составляющая тело среда изотропна, недисперсна и линейна, то, как мы уже показали, задание материального тела заключается в принятии коэффициентами (1.17), (1.18) соответствующих значений в тех узлах сетки, которые содержатся внутри тела.

Источник электромагнитной волны. Самый простой способ задания источника заключается в задании временной зависимости величины Jsource в (1.16). Такой тип источника обычно используется при моделировании диполей. Для генерации плоской волны более удобен другой тип источника, реализуемый с помощью метода полного и рассеянного поля (total-eld / scattered-eld method) [6, 17].

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

Детекторы могут размещаться в произвольном месте вычислительного объема и не привязаны к сетке. Значения полей на них получаются путем интерполяции по соседним сеточным узлам.

Обычный сценарий численного эксперимента FDTD выглядит так:

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

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

• Источник генерирует конечную во времени электромагнитную волну, спектральный состав которой должен покрывать интересующий нас диапазон частот. Далее, волна падает на тела, перерассеивается на них, и через какое-то время уходит в поглощающие граничные условия. История распространения волны фиксируется детекторами.

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

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

Метод полного и рассеянного поля используется для моделирования бесконечно удаленного источника плоской волны. Этот метод основан на линейности уравнений Максвелла и вытекающего из них принципа суперпозиции. А именно, мы предполагаем, что измеряемое (полное) поле Etotal и Htotal может быть представлено в виде суммы Einc, Hinc соотвествует полю падающей волны, которая предполагается известной во всех точках пространства в любой момент времени. Это та волна, которая распространялась бы в пространстве, если бы в нем не существовало никаких тел. Escat, Hscat соответствует рассеянной волне, представляющей собой результат взаимодействия падающей волны с телами.

Значение рассеянной волны заранее неизвестно.

Разностные уравнения FDTD могут независимо применяться как для полного поля, так и для падающего или рассеянного полей. Как показано на Рис. 1.3, это позволяет разбить вычислительный объем на две области: область полного поля, в которой рассчитывается полное поле, и область рассеянного поля, в которой рассчитывается рассеянное поле. Эти две области разделяет виртуальная граница, которая служит для генерации плоской волны в область полного поля. Разностные уравнения, используемые для расчета компонент поля в прилегающих к этой границе сеточных узлах, отличаются от исходных наличием дополнительных слагаемых, в которых учитывается значение поля падающей волны.

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

вычислительного объема тонкого слоя специального материала, называемого идеально согласованным слоем (Perfectly Matched Layer - PML) [6].

Этот материал в идеале полностью поглощает все падающие на него волны без какого-либо отражения независимо от угла падения и длины волны.

Понятие идеально согласованного слоя (PML) было введено Дж.-П.

Беренгером в 1994 году [18].

Работа такого слоя основывалась на разбиении исходных полей E и H на две компоненты, для каждой из которых должны решаться свои уравнения. В последствие были предложены усовершенствованные формулировки PML эквивалентные первоначальной формулировке Беренгера. Так, в одноосном PML (Uniaxial PML) [19, 20] используется анизотропный поглощающий материал, что позволяет не вводить дополнительные переменные и остаться в рамках исходных уравнений Максвелла.

Однако одноосный PML, как и PML в формулировке Беренгера, не удобны тем, что в них отсутствует поглощение затухающих волн, что не позволяет помещать PML близко к рассеивающим телам. Этого недостатка лишен оборотный PML (Convolutional PML) [21, 22], основанный на аналитическом продолжении уравнений Максвелла в комплексную плоскость таким образом, что их решение экспоненциально затухает. CPML также удобнее в ограничении бесконечных проводящих и дисперсных сред. Помимо этого математическая формулировка CPML обладает большей наглядностью и Рис. 1.4. Моделирование с помощью FDTD рассеяния плоской волны на шаре.

доступностью для понимания.

В качестве иллюстрации мы хотим показать ход численного моделирования рассеяния плоской волны на шаре, в котором используются метод полного и рассеянного поля и поглощающие граничные условия (Рис. 1.4). Этот эксперимент интересен тем, что он позволяет сравнить результаты FDTD с известным аналитическим решением этой задачи (решение Ми [23]). Такое сравнение будет приведено в параграфе 1.2.

FDTD может применяться не только для конечных (как на Рис. 1.4), но и для бесконечных периодических структур. Для их моделирования используются периодические граничные условия по одному или по нескольким направлениям. В частности, фотонно-кристаллические пластинки и антиотражающие покрытия, рассматриваемые в главах 3,4, являются планарными периодическими структурами: они периодичны по двум направлениям и имеют ограниченную протяженность по оставшемуся направлению.

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

Рис. 1.5. Схема численного эксперимента FDTD для расчета спектров прохождения и отражения планарной периодической структуры в случае нормально падения. 1 – генерирующая электромагнитный импульс граница между областями полного (справа от 1) и рассеянного (слева от 1) поля; 2, 2’ - плоскости, в которых расположены детекторы для измерения отраженной и прошедшей волны соответственно.

В практических приложениях, как правило, наиболее интересны спектры прохождения и отражения от рассматриваемых структур. Для их получения в случае нормального падения в FDTD используется следующий численный эксперимент (Рис. 1.5). Граница между областями полного и рассеянного поля генерирует плоскую волну в форме ограниченного по времени импульса, который падает на структуру (на Рис. 1.5 этот импульс генерируется на границе 1 и движется слева направо). Часть падающего импульса отражается от структуры, пересекает плоскость, в которой находятся детекторы, меряющие отраженную волну (плоскость 2 на Рис.1.5), и поглощается задним слоем PML. Другая часть проходит сквозь структуру, пересекает плоскость, в которой расположены детекторы, меряющие прошедшую волну (плоскость 2’ на Рис. 1.5), и поглощается передним слоем PML. Численный эксперимент продолжается до тех пор, пока сигнал не выйдет из вычислительного объема. Время выхода сигнала зависит от конкретной геометрии эксперимента. При увеличении продольной протяженности структуры это время увеличивается по причине увеличения расстояния, которое должен пройти сигнал, чтобы из нее выйти. Время выхода сигнала зависит также и от оптических свойств исследуемой структуры, например, при наличии поглощения оно меньше.

По окончании эксперимента записанные детекторами значения полей E(t), H(t) переводятся в частотное представление с помощью дискретного преобразования Фурье. Далее, рассчитывая среднее значение вектора Пойнтинга по времени S() на каждом из детекторов по формуле [23] мы можем посчитать средний поток энергии W (), проходящий через поверхность A, на которой расположены детекторы:

Нормируя поток энергии прошедшей и отраженной волны на падающий импульс, мы получаем интересующие нас спектры прохождения T () и отражения R(). Поглощение вычисляется как A() = 1 T () R().

Описанию численного эксперимента в случае наклонного падения будет посвящен параграф 1.3.

1.2. Метод подсеточного сглаживания Как и в любом другом разностном методе, в FDTD существует проблема неточного отображения границы тела на вычислительную сетку. В непосредственной близости от границы двух сред уравнения Максвелла должны решаться с учетом граничных условий для векторов E и H. Любая кривая поверхность, разделяющая соседние среды и геометрически не согласованная с сеткой, будет искажаться эффектом лестничного приближения“. В результате порядок точности FDTD понижается с изначального второго до первого.

Для решения данной проблемы были предложены различные методы. Первая группа методов основана на изменении способа дискретизации уравнений Максвелла с той целью, чтобы сетка максимально соответствовала рассматриваемой геометрии. Например один из таких методов заключается в увеличении разрешения сетки в тех областях пространства, где расположены тела со сложной геометрической структурой [24]. Также можно видоизменять разностные уравнения в узлах сетки, находящихся вблизи границы между соседними телами [25]. Более радикальным шагом является введение нерегулярной неортогональной сетки, которая полностью согласовывалась бы с рассматриваемой геометрией [26]. Такая сетка автоматически генерируется программой, которая следит за расстановкой всех тел в пространстве. Общим недостатком упомянутых методов является увеличение объема требуемой памяти и порой существенное снижение производительности.

Другая группа методов основывается на введении эффективной диэлектрической проницаемости вблизи границы между телами. Рассмотрим некоторый контрольный объем x y z, окружающий выбранный узел сетки, и предположим, что он содержит границу раздела между двумя средами со значениями диэлектрической проницаемости 1 и 2. Тогда граничные условия для векторов E и D могут быть реализованы на сетке путем ввода тензора обратной диэлектрической проницаемости в форме [7, 27–29]:

где P - это матрица Pij = ni nj соответствующая проектору на вектор нормали n к границе между двумя средами, а есть усреднение по объему.

Выражение (1.24) в рамках FDTD было впервые получено в работе [27] для двумерного случая, но в ней использовались только диагональные элементы матрицы P.

Различные методы введения эффективной диэлектрической проницаемости могут рассматриваться как аппроксимации (1.24). Наиболее простым является применение усредненного по объему или по контуру значения < > [30]. Другие аппроксимации (такие как формула Канеды [31] или VP-EP формула [27]) учитывают направление вектора n. Однако любое частичное использование (1.24) не повышает порядок точности по сравнению с лестничным приближением. Точное выражение (1.24) было использовано в [7, 28, 29], в частности в [7] было показано, что применение тензора (1.24) приводит ко второму порядку точности.

До сих пор метод введения эффективной диэлектрической проницаемости применялся только к диэлектрическим средам. Мы развили схему, предложенную в [7,28,29] также и для случая проводящих и дисперсивных сред.

Рассмотрим произвольный объем, внутри которого пролегает граница между двумя различными средами с зависящей от частоты диэлектрической проницаемостью 1,2 (). Выпишем закон Ампера в частотном представлении, используя выражение (1.24):

где f1 есть доля объема, занимаемая средой с диэлектрической проницаемостью 1 (). Чтобы решить (1.25) во временном представлении, представим E как сумму следующих трех вспомогательных векторов:

Полученные уравнения для E m могут быть решены во временном представлении традиционным для FDTD путем. Например, если зависимость диэлектрической проницаемости от частоты представляется в виде суммы произвольного числа членов Друде или Лоренца, эти уравнения могут быть решены методом вспомогательных дифференциальных уравнений [32]. Если зависимость диэлектрической проницаемости от частоты дана в виде экспериментальной таблицы, то ее можно аппроксимировать необходимым числом членов Друде или Лоренца.

Рис. 1.6. Геометрия численного эксперимента FDTD (масштаб не сохранен).

Мы реализовали описанный метод в качестве составной части написанной нами библиотеки EMTL, о которой пойдет речь в главе 2. Вычисление вектора нормали n к границе между двумя средами и расчет доли объема f1 происходят на начальной фазе работы во время рассмотрения тел, пересекающих контрольный объем, содержащий в себе центр обновляемого контура. Этот контрольный объем обычно представляет собой параллелепипед, равный по объему ячейке Йи. Вспомогательные вектора E m всегда перпендикулярны нормали к плоскости сеточного контура, вследствие чего для каждого такого контура требуется только одна скалярная проекция E m. Диагональные компоненты Pkj ( H)j, j = k вычисляются во время общей процедуры обновления значений в узлах сетки, в то время как недиагональные компоненты j = k получаются путем линейной интерполяции значений поля H по смежным узлам сетки. Добавление недиагональных слагаемых, временной шаг для вспомогательных переменных (1.27)-(1.29) и их сложение (1.26) происходит после общей процедуры обновления сетки, образуя таким образом специальную сглаживающую поправку“.

Для проверки работы изложенного метода мы сравнили полученные с помощью него результаты для рассеяния плоской волны на сфере с результатами, полученными с помощью теории Ми [23, 33]. Геометрия соответствующего численного эксперимента изображена на (Рис. 1.6). Расчет проводился на сетке Йи. Центр сферы находился в центре ячейки Йи.

Для генерации падающего на сферу волнового импульса использовался efficiency factor for scattering Рис. 1.7. Зависимость фактора эффективности рассеяния Qsca от радиуса проводящей сферы с диэлектрической проницаемостью = 40 и проводимостью = 2 на длине волны = 25.

метод общего и рассеянного поля. На границе сетки были заданы поглощающие граничные условия в форме одноосного идеально согласованного слоя (Uniaxial Perfectly Matched Layer, UPML). Вокруг сферы размещались ближние“ детекторы, которые образовывали замкнутую поверхность. Для получения значений поля на виртуальных дальних“ детекторах, расположенных на значительном расстоянии от сферы, применялся метод преобразования ближних полей в дальние [6, 34]. Значения полей на дальних детекторах использовались для сравнения с решением Ми.

Мы проводили тестирование следующих методов:

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

• усреднение диэлектрической проницаемости по контуру или по объему;

• усреднение обратной диэлектрической проницаемости, которое соответствует P = 1 в (1.25);

• VP-EP метод [27], который был нами реализован для проводящих и дисперсивных сред;

Рис. 1.8. Зависимость относительной погрешности расчета Qsca от разрешения сетки. Условия эксперимента те же, что и для Рис. 1.7.

• собственно применение полного выражения (1.24).

В первой серии тестов мы рассматривали рассеяние плоской волны на проводящей сфере ( = 40, = 2) при различных значениях ее радиуса (Рис. 1.7). Шаг сетки по пространству равнялся r = 1. Для сравнения с теорией Ми мы выбрали фактор эффективности рассеяния, который определяется как сечение рассеяния, нормированное на площадь проекции рассеивателя на плоскость, расположенную перпендикулярно направлению падения (для сферы фактор эффективности рассеяния есть Qsca = S/(4r2 ), где S - сечение рассеяния).

Ступенчатый характер полученной зависимости Qsca при разных радиусах сферы для лестничного приближения и усреднения по контуру обусловлен ступенчатым характером изменения коэффициентов в разностных уравнениях в случае пересечения сферой нового сеточного контура при увеличении ее радиуса. Для усреднения по объему полученная кривая является более гладкой. Отсюда можно сделать вывод, что усреднение по объему важно для распознавания очень маленьких объектов, но не приводит к увеличению точности для объектов большего размера. Заметим, что точки, соответствующие результатам для прямого и обратного усреднения Рис. 1.9. Зависимость фактора эффективности рассеяния Qsca от радиуса для свинцовой сферы, = 1.5 µm.

< >1 и < 1 >, расположены по разные стороны от теоретической кривой. Таким образом, применение тензора обратной диэлектрической проницаемости представляет оптимальную комбинацию обоих способов усреднения (на Рис. 1.7 результаты для метода VP-EP и собственно тензорного метода не различимы от теоретических). Сравнение точности разных методов при различном разрешении сетки 1/r представлено на Рис. 1.8, откуда можно видеть, что использование тензора обратной диэлектрической проницаемости приводит к увеличению порядка точности.

Во второй серии тестов мы рассматривали рассеяние плоской волны = 1.5µm на свинцовой сфере, рассчитанное на сетке с шагом r = 50 nm (Рис. 1.9, 1.10). Для диэлектрической проницаемости свинца мы чениями параметров взятыми из [35]. Отметим, что на рассматриваемой длине волны действительная часть диэлектрической проницаемости свинца отрицательна: 81 + 18i. Поскольку решение уравнений Максвелла существенно различается для областей с положительной и отрицательной диэлектрической проницаемостью, методы, не учитывающие направление границы раздела между соседними средами и основанные на частичной апtheory Рис. 1.10. Угловая зависимость элемента матрицы рассеяния S11 [23] для свинцовой сферы радиуса r = 250 nm, = 1.5 µm.

проксимации выражения (1.24), показывают существенное снижение точности по сравнению методом, в котором примененяется тензор обратной диэлектрической проницаемости.

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

1.3. Итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру Одной из задач, для которой может успешно применяться метод FDTD, является моделирование оптических свойств планарных периодических структур. Такие структуры периодичны по двум направлениям и имеют ограниченную протяженность по оставшемуся третьему направлению (фотонно-кристаллические пластинки, массивы антенн и т. д.). Для получения спектров прохождения и отражения таких структур в методе FDTD генерируется падающая на структуру плоская волна в форме ограниченного по времени импульса, которая распадается после взаимодействия со структурой на прошедшую и отраженную волны.

Поле в рассматриваемой волне F (x, t) в любой точке x и в любой момент времени t должно удовлетворять следующему условию:

где u - единичный вектор, определяющий направление падения плоской волны, а a = m1 a1 + m2 a2, m1,2 Z - комбинация векторов примитивных трансляций двумерной решетки a1,2, которые задают периодичность структуры по двум направлениям.

В случае нормального падения направление, задаваемое вектором u, перпендикулярно плоскости периодичности структуры (ua = 0). Вычислительный объем FDTD в этом случае может быть ограничен одной элементарной ячейкой, на границах которой применяются периодические граничные условия означающие, что в выходящих за границу вычислительного объема точках xb значения поля полагаются равными значениям в точках у противоположной относительного периодического направления границы xb + ab (ab = ±a1,2 для четырех периодических границ).

Если волна по отношению к плоскости периодичности структуры падает наклонно, в граничные условия должен быть включен временной сдвиг где 1,2 есть величина угла падения в плоскости, образованной соответствующим направлением трансляции двумерной решетки и нормалью к ней.

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

Для решения этой проблемы в литературе было предложено несколько различных методов, которые можно условно разделить на три группы. В первой группе вводятся новые переменные, выражаемые через поля E и H таким образом, чтобы устранить временную задержку между элементарными ячейками. В новых переменных уравнения для всех элементарных ячеек эквивалентны, и к ним применимы периодические граничные условия. Однако эти уравнения отличны от уравнений Максвелла, и для их решения требуется написание новых алгоритмов, отличных, например, от алгоритма Йи. Наиболее удачными алгоритмами дискретизации этих уравнений оказались алгоритм введения дополнительных сеток [36] и алгоритм разделения полей [37]. В обоих этих алгоритмах вводятся дополнительные переменные, что приводит к увеличению объема требуемой памяти.

Более весомым их недостатком является увеличение времени расчета, которое обусловлено уменьшением максимального для устойчивости метода значения t/x при увеличении угла падения. Для углов близких к требуемое максимальное отношение t/x стремится к нулю, а расчет становится невозможным.

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

в которых временной сдвиг заменяется приращением фазы b. Такая замена допустима в случае падающей монохроматической волны частоты (как это происходит в методе синусов и косинусов [38]), где фаза связана с углом(ами) падения с помощью соотношения b = ±0 a1,2 sin 1,2 /c.

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

Однако, в качестве падающей волны можно использовать пакет плоских волн, падающих под разными углами, но имеющих общее значение фазы b, как это предлагается в методе спектрального FDTD [39]. Результаты такого эксперимента будут относится уже к диапазону пар (, ). Недостатком описанных двух методов является то, что они фактически переходят в частотное представление, что лишает их тех преимуществ FDTD, которые являются следствием его работы во временном представлении.

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

В методе углового цикла [6] также используются дополнительные ячейки. Сам метод основан на неодновременном обновлении значений полей на сетке, что позволяет получать эти значения из будущего, пользуясь значениями в соседних опережающих по времени ячейках. Недостаток этого метода в том, что он применим только для относительно небольших углов (в трехмерном случае допустимы только углы меньшие 350 ).

В этой главе мы предлагаем новый итерационный метод для моделирования наклонного падения плоской волны на периодическую планарную струкутру. Этот метод основан на итерационном приближении к искомому решению, удовлетворяющему условиям (1.30). Для моделирования падения плоской волны под выбранным углом требуется несколько (5-10) отдельных экспериментов, которые мы в дальнейшем будем называть итерациями, и дополнительная память для хранения данных, которые передаются от предыдущей итерации к последующей.

Геометрия численного эксперимента представлена на Рис. 1.11 (в целях простоты изложения мы рассматриваем двумерный случай с примитивными векторами трансляции решетки a1 = a, a2 = 0 и временным Рис. 1.11. Геометрия численного эксперимента. 1, 2 – генерирующие (TF/SF) границы, для которых в качестве источника берется сигнал, переносимый с точек образов 1, 2 ; 3 – TF/SF граница для генерации падающей плоской волны; 4 – детекторы для записи прошедшего сигнала.

сдвигом t = |tb | = a sin /c). Мы используем метод общего и рассеянного поля (Total Field / Scattered Field Method, TF/SF) для генерации плоской волны в область общего поля (эта область заштрихована на Рис. 1.11).

Как полное, так и рассеянное поле в искомом решении удовлетворяет условию (1.30). Зависящие от времени граничные условия TF/SF на границе 3 соответствуют падающей наклонно плоской волне и известны аналитически. Основная идея нового метода заключается в применении дополнительных граничных условий на границах 1 и 2, которые работают подобно границе TF/SF, где в качестве генерируемой волны используются значения полей в точках-образах на границах 1 и 2 (x1 = x1 + a и x2 = x2 a) с соответствующим временным сдвигом. Для отрицательного временного сдвига в прошлое могут быть использованы значения полей с текущей i-ой итерации, а для положительного сдвига в будущее мы используем значения полей, записанные на предыдущей i 1-ойитерации:

Рис. 1.12. Распределение электромагнитной энергии на первой и пятой итерации численного эксперимента с наклонным падением плоской волны на металлическую пластинку, = 450.

Падающая и отраженная волны отчетливо видны на пятой итерации.

Расстояние между границами 1 и 2 задается больше периода a на некоторое расстояние x, равное нескольким сеточным шагам, для того чтобы разделить точки-образы и дополнительную границу TF/SF. Для записи значений E(x1, t) и H(x1, t) на границе 1 используется дополнительный буфер памяти размером 6N T /t, где T - продолжительность одной итерации, t - шаг по времени, а N = Ny Nz - количество сеточных узлов, приходящихся на границу 1 (N T /t для каждой из трех компонент полей E и H). Для границы 2 требуется буфер размера 6N t/t. Для t < 0 и t > T поля в точках-образах полагаются равными нулю Fi (x1,2, t) = 0. На первой итерации в качестве граничного условия на границах 1 и 2 берется аналитическое значение для падающей плоской волны или нулевой сигнал.

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

На Рис. 1.12 продемонстрирована сходимость итерационного метода к искомому решению. Как видно из рисунка, исчезновение рассеянного поля слева от границы 1 и справа от границы 2 свидетельствует о получении искомого численного решения. Для измерения отличия значения волны на текущей итерации от искомого решения может использоваться величина потока энергии в область рассеянного поля через границы 1 и 2 (Рис.1.13).

Для доказательства сходимости метода представим генерируемый на i-ой итерации на границах 1 и 2 сигнал в виде суммы Fi = F0 + Fi, где Рис. 1.13. Поток энергии через границы 1 и 2, соответствующий численной ошибке. Эксперимент соответствует Рис. 1.12.

F0 (xb, t) удовлетворяет условию (1.30), а Fi есть сигнал-ошибка. По причине линейности уравнений Максвелла и операции временного сдвига по отношению к полям Fi (x, t), а также вытекающего из линейности принципа суперпозиции, эволюция ошибки Fi при увеличении i может рассматриваться отдельно от решения F0 (xb, t). На каждой следующей итерации сигнал-ошибка генерируется на периодической TF / SF границе с (а) временным сдвигом в будующее и (б) меньшей амплитудой, как это видно на Рис.1.13. Для непосредственного доказательства этого утверждения рассмотрим некоторый сигнал Fi, генерируемый границей 1 в течение некоторого интервала времени, начинающегося в t = 0. Тогда:

• Время, необходимое для того, чтобы сигнал распространился до точки-образа 1 не меньше t11 = |x1 x1 |/c = a/c. В результате этого сигнал-ошибка будет записан на 1 и перенесен на границу 1 с положительным временным сдвигом t согласно соотношению (1.35), то есть общее время задержки равно t1 = t11 t = a(1 sin )/c, что для углов = 900 больше нуля. Рассматриваемый сигнал также будет записан на границу 2 и возникнет на границе 2 на следующей итерации, но будет распространяться в противоположную сторону от направления к границам 1 и 1. На следующей итерации, по крайней мере, в течение временного интервала [0, t1 ) на границе 1 сигналошибка Fi+1 будет равен нулю. Эти же рассуждения применимы к границе 2, при этом временная задержка t2 = t1 + 2x/c. Таким образом, любой генерируемый сигнал на границах 1 или 2 временной длительности T и содержащий в себе ошибку, не удовлетворяющую условию (1.30), будет отделен от этой ошибки по крайней мере за n T /t1 итераций (за это число итераций ошибка запаздает на время T ).

В приведенных рассуждениях мы предполагали, что аналитическая граница 3 не генерирует численной ошибки, что в общем случае не верно. На Рис.1.13 видны два ошибочных сигнала. Первый из них это небольшой стабильный во времени сигнал, который сопровождает падающую волну и связан с численной ошибкой метода TF/SF.

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

Так, на Рис.1.13 уже на 20-ой итерации это сигналы оказываются полностью разделены по времени.

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

Мы реализовали описанный метод в качестве составной части библиотеки EMTL, о которой пойдет речь в главе 2. Метод был протестирован при Рис. 1.14. Зависимость коэффициента прохождения однослойного фотонного кристалла, состоящего из упакованных в квадратную решетку металлических шариков ( = 40, = 2, период решетки a = 1, радиус сфер 0.375) от частоты падающей волны, = 450.

моделировании однослойного фотонного кристалла, состоящего из упакованных в квадратную решетку металлических шариков ( = 40, = 2), на который падал ограниченный во времени волновой импульс с плоским фронтом. Период квадратной решетки полагался равным единице a = 1, радиус сферы 0.375, а шаг по сетке 0.05 (c = 0 = 1). Для уменьшения численной ошибки, связанной c неточным отображением границ шариков на сетку, использовался метод подсеточного сглаживания [7, 8].

На Рис. 1.14, 1.15 представлено сравнение результатов для спектра прохождения, полученных с помощью описанного метода и с помощью метода LKKR [41]. Результаты FDTD находятся в хорошем согласии с результатами LKKR уже после пяти итераций, для больших частот и малых углов падения результаты сходятся быстрее.

В заключение мы хотим отметить практические отличия нового метода от ранее предложенных методов. Новый метод лишен таких недостатков, как неустойчивость при отвесных углах падения, увеличение требуемого размера сетки, переход в частотное представление. С помощью нового метода может моделироваться падение плоской волны любой длительности по времени под произвольным углом, однако иногда это может потребовать увеличения требуемого объема памяти и числа итераций (а, следовательно, Рис. 1.15. Зависимость коэффициента прохождения однослойного фотонного кристала (параметры фотонного кристалла те же, что и для Рис. 1.14) от угла падения, f = 0.6/a.

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

1.4. Выводы В рамках метода FDTD предложен и реализован алгоритм подсеточного сглаживания для дисперсных сред, позволяющий существенно уменьшить величину численной ошибки, связанной с неточным отображением тел, обладающих произвольной формой поверхности, на прямоугольную вычислительную сетку. Точность метода проанализирована путем сравнения полученных с его помощью результатов с аналитическим решением Ми задачи о рассеянии плоской волны на сфере. Показано, что его применение приводит к увеличению порядка точности метода FDTD по сравнению с обычным лестничным“ приближением.

Предложен и реализован итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру в численном расчете FDTD. С его помощью можно получать результаты во всем диапазоне углов от 00 до 900, а также наблюдать реальное поведение поля во времени. Особенностью метода является то, что получение решения для одного угла требует нескольких отдельных итераций-экспериментов. Точность метода проанализирована путем сравнения получаемых с помощью него результатов с результатами, полученными с помощью слоевого метода Коринга-Кона-Ростокера.

2. Применение шаблонного метапрограммирования для эффективной реализация FDTD В первоначальном узком смысле под FDTD подразумевалось использование базового алгоритма Йи для численного решения уравнений Максвелла. В современном более широком смысле FDTD включает в себя множество самых разнообразных возможностей: моделирование сред с дисперсными и нелинейными свойствами, применение различных типов сеток (помимо первично предложенной прямоугольной сетки Йи), использование методов постпроцессорной обработки результатов и т. д. Для реализации всех этих возможностей требуется написание больших параллельных программ, которые обладали бы высокой производительностью, были бы легко расширяемы и обладали доступным, читаемым кодом. К сожалению, распространённые некоммерческие программы FDTD (например, MEEP, XFDTD, CFDTD) лишены этих свойств и не предназначены для широкого круга задач. Это явилось причиной разработки собственной модели FDTD, которую мы будем в дальнейшем называть EMTL (Electromagnetic Template Library).

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

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

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

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

EMTL реализована в виде библиотеки шаблонов C++, которые могут быть использованы в качестве строительных компонент при создании конкретных численных приложений. Для реализации возможности проведения параллельных расчетов используется библиотека MPI (Message Passing Interface) [43].

Библиотека EMTL находится в открытом доступе в интернете (http://fdtd.kintechlab.com).

2.1. Контурный подход к дискретизации уравнений Максвелла В качестве отправной точки для описания упомянутого в названии параграфа контурного подхода служит первая пара уравнений Максвелла в интегральной форме, которая в отсутствие зарядов и в случае изотропных, недисперсных и линейных сред (см. (1.5), (1.6)) имеет вид Рис. 2.1. H-контур, служащий для обновления поля E.

Рассмотрим некоторый контур, ограниченный выпуклым многоугольником (Рис.2.1). Аппроксимируем поверхностный и контурный интегралы в (2.1) и (2.2) следующим образом:

где для электрических и магнитных постоянных используются обозначения perm1 =, loss1 =, perm2 = µ, loss2 =. Поля F1,2 = E, H в (2.3) находятся по разные стороны от знака равенства: поле F c и ток J c, измерямые в центре контура и умножаемые на его площадь S, находятся слева, а поле F l, измеряемое в центральных точках сторон контура и умножаемое ветора l, находится справа.

Как это делалось в параграфе 1.1.1, используя для поля в левой части (2.3) центральное разностное отношение для производной по времени и линейную интерполяцию по соседним временным слоям и явно выражая в (2.3) значение поля F c на следующем временном шаге, получаем:

F1,2 (t + t)S = 1,2 F1,2 (t t)S + 1, где 1,2, 1,2, 1,2 - некоторые коэффициенты (для сетки Йи эти коэффициенты задаются выражениями (1.17), (1.18)).

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

В случае, когда контур пересекает границу между двумя различными средами (Рис. 2.1), для корректного задания электромагнитных свойств в левой части уравнения (2.3) необходимо использовать какое-либо приближение. Самое простое приближение соответствует заданию этих свойств в центре контура ("лестничное приближение"). Более точным является применение тензора обратной диэлектрической проницаемости, о котором речь шла речь в параграфе 1.2.

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

где I(x) есть множество индексов (обычно небольшое), определяемое точкой x, а ci - интерполяционные коэффициенты (Рис.2.2). Операция Li связывает i-ое слагаемые (i-ая точка на Рис.2.2) с массивом данных. Обычно эта операция ставит в соответствие каждой точке уникальную ячейку памяти: Li (M ) = Mi. Возможны более изысканные случаи, например, когда Рис. 2.2. Интерполяций значений поля на опорных точках контура по узлам сетки.

значения в точке получаются путем вычисления интерполяционных полиномов.

В качестве примера рассмотрим, как осуществляется линейная интерполяция на трехмерной ортогональной сетке N0 N1 N2 с шагом по пространству x. В этом случае значение для интерполяции проекции поля на выбранное направление в произвольной точке используются слагаемых: 8 вершин содержащего эту точку куба, задаваемого соседними узлами сетки, умножить на 3 компоненты поля для каждого узла. Если 3 компоненты поля нумеровать индексом k, а вершины куба нумеровать тройкой индексов m0, m1, m2, принимающих значения либо 0, либо 1, то тогда в случае, когда данные для каждой компоненты хранятся в памяти вместе, множество индексов из (2.7) есть i8k+4m0 +2m1 +m2 = N0 N1 N2 k + N1 N2 (n0 + m0 ) + N2 (n1 + m1 ) + n2 + m2, (2.9) где приращение индекса в массиве данных на единицу соответствует изменению вдоль оси z (j = 2). Более подробно о способах хранения данных в памяти будет идти речь в параграфе 2.4.1.

Каждой компоненте ставится в соответствие уникальная ячейка памяти:

Интерполяционные коэффициенты равны где nj = xj mod x суть номера сеточных узлов, расположенных вблизи точки x, а ek - единичные вектора вдоль соответствующих направлений.

Выражение (2.7) может быть использовано как для суммирования входных значений поля перед обновлением выходного значения поля в (2.6), так и для помещения в память выходного значения поля. Пользуясь линейностью (2.6) и (2.7) по отношению к Li (M ), а также предполагая наличие линейности у операции Li, можно показать, что каждая ячейка памяти может обновляться явным образом через значения в остальных ячейках. Это позволяет на основании (2.6), (2.7) реализовать цикл обновления поля в самом общем виде, не зависящем от типа используемой сетки. Однако, с учетом организации конкретной сетки выражение (2.7) может быть сильно упрощено, в результате чего число входных и выходных слагаемых в (2.6) сократится. Например, в сетке Йи множество векторов I(x) всегда равно одному, а для обновления значения поля по контуру можно пользоваться простым выражением (1.16), использующем лишь одну выходную и пять входных ячеек.

2.2. Код как совокупность взаимодействующих компонент EMTL можно представить как совокупность взаимодействующих друг с другом компонент, причем базовые примитивные компоненты могут быть задействованы в качестве параметров шаблонов при построении компонент более высокого уровня. Такими базовыми компонентами являются контура, итераторы (они работают подобно итераторам из стандартной библиотеки шаблонов C++ [44]), множества интерполяционных коэффициентов, сетки и т. д. В целях увеличения производительности методы базовых компонент, постоянно вызываемые в основных внутренних циклах, реализованы в виде inline функций.

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

Соответствующие телам компоненты наследуют от своего базового класса SpaceRegion следующую функциональность:

• умение определять, находится ли определенная точка внутри тела;

• нахождение площади сечения телом заданного контура;

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

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

Ключевым строительным блоком EMTL является сетка или класс emMeshBlock. Сетка задает совокупность контуров, участвующих в разностной схеме, организует размещение данных в памяти и определяет способ интерполяции значений поля по сеточным узлам в своих внутренних точках. Более подробно, сетка должна уметь:

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

• создавать интерполяцию в виде (2.7) для произвольной внутренней • возвращать итератор по всем E- и H-контурам, формирующим сетку;

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

В рамках указанных требований могут быть реализованы сетки самого разного типа, включая не ортогональные и неоднородные. Заметим, что совокупность контуров, формирующих сетку, должна быть согласована со второй парой уравнений Максвелла (1.3,1.4). Поскольку в EMTL это условие не проверяется, за него целиком отвечает создатель сетки.

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

2.3. Стадии расчета В предлагаемой модели FDTD можно выделить следующие три стадии расчета:

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

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

• анализ записанных данных.

2.3.1. Задание геометрии На стадии задания геометрии происходит перебор всех контуров, формирующих сетки, используемые в данном численном эксперименте. На основании анализа расположения этих контуров относительно заданных тел рассчитываются коэффициенты в соответствующих разностных уравнениях. Эта процедура выполняется только один раз перед вычислительным циклом, и время, затрачиваемое на нее, существенно меньше времени, затрачиваемого на сам вычислительный цикл. Таким образом, производительность этой процедуры не столь важна, и ее можно реализовать в наиболее общем виде, не зависящем от внутреннего устройства используемых сеток. Для того, чтобы это было возможно, сетка должна предоставлять итератор по всем своим внутренним контурам, которые в свою очередь должны:

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

• возвращать последовательность всех своих вершин;

• устанавливать соответствие между опорными“ точками контура (центр контура и центральные точки его сторон) и соответствующими коэффициентами и индексами в массивах компоненты emMeshBlock;

• предоставлять некоторый "контрольный объем“, содержащий в себе данный контур.

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

Участвующие в задании геометрии компоненты помещаются в объект типа emBlockContainer, через который осуществляется управление ходом численного эксперимента. Этот объект может включать в себя произвольное число сеток, каждая из которых должна принадлежать одному из двух заранее указанных типов.

Рис. 2.3. Контур, пересекающий границу между двумя сетками.

Контейнер должен содержать ограничивающую область, для которой задаются соответствующие граничные условия. Эти условия определяют алгоритм действий, который применяется, когда определенный контур пересекает ограничивающую область, так что одна или несколько из его входных точек находятся снаружи ее. В случае периодических граничных условий значения полей в этих точках берутся с противоположной стороны ограничивающей области, в случае же нулевых“ граничных условий - полагаются равными нулю. Если расчет ведется для одномерной или двумерной задачи, на дополнительные измерения накладываются периодические граничные условия, а размер вычислительного объема по этим измерениям сокращается до одного - двух контуров.

Каждой сетке ставится в соответствие целое число. Это число определяет приоритет обращения к той или иной сетке при их взаимном пересечении (как правило, сетки с большим разрешением обладают большим приоритетом). А именно, при запросе у контейнера значения поля в определенной точке пространства, если эта точка находится в области пересечения нескольких сеток, значение поля будет интерполироваться в первую очередь по узлам сетки, обладающей более высоким приоритетом. Сетки с более низким приоритетом используют значения полей, получаемые у сеток с более высоким приоритетом, в качестве входных точек своих контуров, которые расположены в области пересечения сеток (Рис. 2.3).

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

Каждый такой алгоритм может быть организован последовательной обработки нескольких отдельных потоков данных. Сетке Йи соответствует минимальное число потоков: 5 для входных данных и 1 для выходных.

Как правило, каждый раз, когда итератор по контурам достигает границы сетки, поток пересвязывается с другим адресом в памяти.

Моделирование дисперсных сред, источников, поглощающих условий и т. д. требует модификации базового алгоритма (2.6). Традиционном подходом является включение специализированных алгоритмов, реализующих всевозможные модификации (2.6), непосредственно в основной цикл. Каждому такому алгоритму ставится в соответствие одно из возможных значений некоторого флага. В основном вычислительном цикле этот флаг используются для определения того, какой именно алгоритм использовать для каждого контура.

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

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

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

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

Как уже упоминалось в параграфе 1.1.1, максимально возможная величина шага по времени t в FDTD связана с шагом по пространству x и скоростью света c условием Куранта: t x/c. По этой причине сетки различным пространственным разрешением должны обновляться с различным временным шагом.

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

2.4. Оптимизация использования памяти Последовательность контуров, к которым применяется та или иная поправка, определяется рассматриваемой геометрией и обычно содержит в себе некоторую закономерность. Помимо этого часто используемые поправками к разным контурам данные (например, коэффициенты в дополнительных разностных уравнениях) совпадают. Поэтому в целях экономии оперативной памяти в EMTL используется простой способ упаковки этих данных. Он основан на выделении линейно возрастающих подпоследовательностей в последовательности целых чисел. Каждая такая подпоследовательность заменяется на тройку чисел: целый код -1, величина приращения и число элементов в подпоследовательности. Этот способ упаковки достаточно прост, эффективен и не требует больших затрат при распаковке.

Указанным способом могут кодироваться как последовательности итераторов контуров, так и последовательности индексов какого-либо массива данных (например, массива коэффициентов в разностных уравнениях).

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

2.4.1. Тестирование различных способов организации в памяти Во всех разностных методах основной является процедура обновления массива сеточных данных. Для организации этого массива может использоваться принцип близости точек в пространстве: данные, соответствующие физическим свойствам в данной точке в пространстве, хранятся вмеPowerPC, objects Рис. 2.4. Время, затрачиваемое на обновление одной ячейки Йи, для различных способов хранения сеточных данных в памяти и для разных размерностей сетки по трем направлениям (см. комментарий в тексте) в случае вакуумной задачи маленького размера (размер сетки MB). Код был скомпилирован на архитектуре PowerPC с помощью компилятора IBM xlc 7. (-q64 -O5) и на архитектуре Itanium с помощью компилятора Intel icc 9.0 (-O3). Тесты проводились на кластере Межведомственного суперкомьютерного центра Российской Академии Наук (МСЦ РАН).

Рис. 2.5. Результаты тестов аналогичным тем, что представлены на Рис. 2.5, для задачи большого размера (размер сетки 1.3 GB) сте.

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

Поскольу в EMTL проведение расчета с использованием сетки и организация внутреннего устройства самой сетки отделены друг от друга, это позволяет отдельно исследовать эффективность двух указанных выше способов организации данных в памяти ("objects"и "pencils"на Рис. 2.4, 2.5).

Для наших тестов были выбраны две ортогональные сетки Йи размером nx ny nz: в одной сетке в смежных ячейках памяти хранилось по 6 компонент полей E и H, относящихся в одной ячейке Йи (Рис. 1.1), а в другой данные для каждой компонеты хранились вместе согласно (2.9). Для обоих случаев наиболее быстро меняющийся пространственный индекс соответствовал направлению z. Вариируя значение nz при постоянном размере сетки и соблюдении условия nx = ny, мы измеряли время, затрачиваемое на обновление одной ячейки Йи. Для обеих архитектур (Itanium 2 и PowerPC 970), на которых производились тесты, второй способ хранения данных оказался эффективнее. Отметим, что величина nz/nx, определяющая количество пересвязываний обрабатываемого потока данных, также влияет на производительность, но в существенно меньшей степени.

2.5. Проведение параллельных расчетов В EMTL поддерживается возможность проведения параллельных расчетов посредством разбиения вычислительного объема на домены. В качестве доменов могут использоваться произвольные объекты класса SpaceRegion. На стадии перебора контуров выявляются граничные контура, входные точки которых принадлежат не содержащим их домену или сетке. В "пересылочные поправки"для таких контуров записываются интерполяционные коэффициенты и индексы в соответствующих массивах данных, с помощью которых на каждом временном шаге определяется значение поля во входных точках таких контуров. Единственным отличием параллельных пересылок между доменами от пересылок между различными сетками внутри одного домена является то, что в первом случае данные пересылаются путем посылки MPI сообщения, а во втором используется прямое копирование данных в оперативной памяти.

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

Разбивающая плоскость перпендикулярна направлению, вдоль которого протяженность объема максимальна. Размеры образуемых в результате бисекции объемов выбираются таким образом, чтобы каждому их них соответствовало одинаковое значение w(x)dV /ni, где ni - число процессоров, приходящихся на i-ый объем, w(x) - "трудозатратность"на единицу объема, которая может быть определена экспериментально для каждого типа Рис. 2.6. Зависимость времени выполнения программы от числа процессоров в серии масштабируемых тестов на кластере ФМБФ МФТИ (Athlon MP 2400, Gigabit). Начальный размер задачи (сетка 1403 ) масштабируется в одном измерении пропорционально числу процессоров.

Количество временных итераций равно 100.

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

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

Например, в случае прямоугольной сетки, когда соответствие между элементами сетки и индексами в массиве данных задается согласно (2.9), протяженность доменов вдоль направления z, которое соответствует наиболее быстро изменяющемуся индексу в массиве данных, выбирается максимально возможной, в то время как разбиение осуществляется по оставшимся двум измерениям.

UNSCALED

SCALED

Рис. 2.7. Результаты для серии параллельных тестов на кластере МСЦ РАН (PowerPC 970, Myrinet). Базовый размер задачи равен 3001002. Эффективность основного вычислительного цикла (без учета параллельных пересылок) изображена залитыми квадратиками.

2.5.1. Измерение параллельной эффективности С целью измерения параллельной эффективности нашего кода мы провели серию масштабируемых и немасштабируемых тестов на двух различных архитектурах (Рис. 2.6, 2.7). В немасштабируемых тестах измеряется зависимость времени работы программы t(N ) от числа процессоров N при сохранении размера сетки постоянным, при этом параллельная эффективность полагается равной t(1)/N t(N ). В масштабируемых тестах размер сетки увеличивается пропорционально числу процессоров, а параллельная эффективность равна t(1)/t(N ).

Для сравнения мы использовали некоммерческую программу MEEP, разработанную в Массачусетском технологическом институте (MIT).

Как видно из Рис. (2.6) время выполнения программы на используемой архитектуре зависит от режима загруженности сети, а это означает, что пропускная способность сетевого коммутатора лимитирует сверху производительность программы. Отметим, что различие во времени выполнения для MEEP и EMTL находится в обратной зависимости от пропускной способности коммутатора, из чего следует, что оно связано главным обРис. 2.8. Ускорение на единицу моделируемого объема для различных архитектур и реализаций MPI, нормировано на время исполнения на 64 ядрах. HT - Harpertown+Inniband, HT - Nehalem+Inniband (измерения в ИНТЕЛ), Xeon - кластер Т-платформы, узлы 2xXeon 2.33 GHz 8 GB RAM, Inniband.

разом со способом распараллеливания. Большой разброс, наблюдаемый в результатах другой серии тестов для типичной задачи на 1000 временных шагов (Рис. 2.7), также свидетельствуют о непосредственной зависимости времени работы программы от сетевой нагрузки системы.

При поддержке компании Intel были проведены тесты на системах Harpertown и Nehalem на большом числе вычислительных ядер (вплоть до тысячи). В этих тестах размер сетки полагался постоянным или зависел от числа процессоров таким образом, чтобы общий объем задачи сохранялся при его увеличении. Задача представляла собой численное решение уравнений Максвелла в диэлектрике и металле. Результаты тестов приводятся на Рис. 2.8, где показано ускорение на единицу моделируемого объема, характеризующее масштабируемость программы. Ускорение вычисляется относительно задачи, исполняемой на 64 процессорах при размере сетки 2003. Оба типа тестов показали линейный или суперлинейный характер масштабируемости приложения. На этом же графике приведены результаты тестов на малом кластере, состоящем из 24 2-процессорных узлов с четырехядерными процессорами Xeon (2xXeon 5345 2.33 GHz 8 GB RAM).

Из этих результатов следует, что характер производительности приложения зависит от архитектуры системы и реализации MPI.

2.6. Выводы Предложен новый подход программной реализации метода решения уравнений Максвелла в конечных разностях (FDTD), основанный на интегральном представлении уравнений Максвелла. В рамках этого подхода написана параллельная программа на языке программирования С++. В ходе ее написания реализован ряд оригинальных решений: оптимизация использования памяти с помощью упаковки используемых в разностных уравнениях коэффициентов; увеличение производительности посредством разбиения алгоритма обновления сеток на основной цикл, поддающийся векторизации, и специализированные поправки к нему; увеличение параллельной эффективности путем балансировки доменов. Специально проведенные тесты показали, что быстродействие программы примерно в два раза выше, чем у программы MEEP, разработанной в Массачусетском технологическом институте. Также был продемонстрирован линейный характер масштабируемости программы вплоть до нескольких тысяч вычислительных ядер. Развитый при написании программы подход обеспечивает возможность расчета с одновременным использованием нескольких сеток, каждая из которых максимально приспособлена к соответствующей ей области моделируемой структуры. Это является особо важным при моделировании оптических свойств структур сложной формы с неоднородностями, много меньшими длины волны. Программа находится в открытом доступе в интернете (http://fdtd.kintechlab.com).

3. Металлические фотонные кристаллы как источники 3.1. Обзор литературы Фотонные кристаллы (ФК) представляют собой структуры, характеризующиеся периодическим изменением диэлектрической проницаемости в пространственных направлениях [2, 45–47]. В диапазоне длин волн, сравнимых с характерным размером изменения диэлектрической проницаемости, оптические свойства ФК сильно отличаются от оптических свойств сплошных сред. Распространение излучения внутри ФК благодаря свойству периодичности становится похожим на движение электрона внутри обычного кристалла под действием периодического потенциала. В результате электромагнитные волны в ФК имеют зонный спектр и координатную зависимость, аналогичную блоховским волнам электронов в обычном кристалле. При определенных условиях в зонной структуре ФК образуются щели, аналогично запрещенным электронным зонам в естественных кристаллах.

Перераспределение плотности фотонных состояний в ФК приводит к изменению спектра их теплового излучения [3]. Внутри запрещенной зоны плотность фотонных состояний равна нулю, и излучение электромагнитных волн подавлено. Вне запрещенной зоны интенсивность этого излучения непосредственно связана с плотностью фотонных состояний и может сильно варьироваться на разных длинах волн. На основании этого в литературе было высказано предположение о возможности практического использования ФК в качестве высокоэффективных источников видимого света [3,48,49]. Излучение в инфракрасном диапазоне (ИК) таких источников должно быть подавлено, что привело бы к увеличению световой отдачи по сравнению с обычными лампами накаливания, примерно 95% излучения которых находится в ИК области, и лишь 5% приходится на видимый диапазон. В связи с этим становится актуальным исследование излучательных характеристик ФК при высоких температурах ( 2400K) и возможности создания на их основе источников света.

В настоящее время имеется лишь ограниченное число работ (см., напр., [48, 49]), в которых исследуется излучательная способность ФК в видимом диапазоне. Вместе с тем, значительная часть результатов, касающихся влияния геометрических параметров ФК на спектр его излучения, получена для ИК диапазона в рамках исследований возможности применения ФК для управления спектром теплового излучения (см. [50–58]). Это отразилось на характере данного литературного обзора. Вначале нашего изложения мы укажем на существующие методы расчета излучательной способности ФК. Затем мы изложим несколько имеющихся работ, в которых исследуется излучение ФК в видимом диапазоне. Далее значительная часть обзора будет посвящена работам, касающихся излучательной способности ФК в ИК диапазоне, поскольку полученные в этих работах результаты имеют общий характер.

Перейдем к способам расчета излучательной способности. Существующие способы можно разделить на две группы: прямые и косвенные [59,60].

Прямое вычисление излучательной способности тела предполагает его разбиение на маленькие излучающие объемы, учет поглощения и перерассеяния излученной каждым объемом энергии на остальных объемах и суммирование по всем объемам. Более подробное описание этого способа можно найти в работе [60], в которой рассчитывалось тепловое излучение нагретой пластинки помещенной на бесконечную подложку. Отметим также работу [61], в которой прямой подход был применен с использованием метода FDTD при исследовании двумерного ФК, состоящего из бесконечных стержней, упакованных в квадратную решетку.

Излучательная способность может быть также рассчитана косвенно с помощью закона Кирхгофа [52,59]. Эквивалентность прямого и косвенного методов была численно продемонстрирована в работах [60, 61]. Поскольку в наших собственных расчетах применяется второй метод, мы опишем его подробнее.

Пусть имеется произвольное излучающее тело. Поток излучаемой им световой энергии, заключенной в частотном интервале d и телесном угле d, в направлении (, ) через единичную площадку перпендикулярную к этому направлению, равен [52] где E(,,, T ) есть излучательная способность тела, а u(, T ) есть спектральная плотность излучения абсолютно чёрного тела [62] Далее, поглощательная способность A(,,, T ) есть доля поглощаемой телом энергии при падении нее излучения, заключенного в частотном интервале d и телесном угле d, в направлении (, ).

Согласно закону Кирхгофа [63–65] В случае плоскопараллельного ФК слоя поглощательная способность равна [52] где T(,,, T ) есть доля прошедшей энергии, а R(,,, T ) - доля отраженной. В используемых нами выражениях подразумевается усреднение по обоим поляризациям света.

Отсюда получаем, что поток излучаемой ФК слоем световой энергии, заключенный в частотном интервале d, через единичную площадку, параллельную плоскости ФК слоя, равен [52] где интегрирование ведется по телесным углам, задающим направление наружу“ из ФК пластинки. Поток энергии во всем частотном диапазоне равен J(, T )d.