Параллельный алгоритм для решения трехмерных уравнений Максвелла с разрывной диэлектрической проницаемостью на тетраэдральных сетках Текст научной статьи по специальности «Математика»

УПРАВЛЕНИЕ, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И ИНФОРМАТИКА

УДК 519.63

Т. З. Исмагилов

ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ РЕШЕНИЯ ТРЕХМЕРНЫХ УРАВНЕНИЙ МАКСВЕЛЛА С РАЗРЫВНОЙ ДИЭЛЕКТРИЧЕСКОЙ ПРОНИЦАЕМОСТЬЮ НА ТЕТРАЭДРАЛЬНЫХ СЕТКАХ

Предлагается конечно-объемный метод для численного решения трехмерных уравнений Максвелла с разрывной диэлектрической проницаемостью на неструктурированных сетках. Метод имеет второй порядок аппроксимации по времени и пространству для разрыва диэлектрической проницаемости, проходящего по произвольной гладкой поверхности. Численный алгоритм допускает параллельную реализацию с помощью геометрической декомпозиции для использования на многопроцессорных ЭВМ. Приведенные результаты тестовых расчетов подтверждают второй порядок предлагаемого метода и высокую эффективность параллельной реализации. Уравнения Максвелла; метод конечных объемов; схема Годунова; разрывная диэлектрическая проницаемость

ВВЕДЕНИЕ

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

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

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

Контактная информация: ismagilov@academ.org Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии-2010».

К настоящему времени было предложено несколько конечно-объемных алгоритмов для решения уравнений Максвелла. В алгоритмах, предложенных в работе [5], электрические и магнитные поля аппроксимируются на смещенных сетках, как и в работе [1]. В алгоритмах, рассмотренных в работах [6-8], все компоненты электромагнитного поля аппроксимируются в центрах ячеек.

С точки зрения приложений, важным является возможность распараллеливания численных алгоритмов, что позволяет применять их на многопроцессорных ЭВМ. Конечно-объемные алгоритмы, использующие неструктурированные сетки, могут быть распараллелены с помощью геометрической декомпозиции. В работе [9] рассматривалась параллельная реализация конечно-объемного алгоритма для решения двумерных уравнений Максвелла.

Одной из основных трудностей при построении схем второго порядка для уравнений Максвелла остается случай разрывных свойств среды. В работе [2] рассматривались различные способы сглаживания разрывной диэлектрической проницаемости. Но ни один из них не позволял сохранять порядок аппроксимации исходной схемы.

Для решения этой проблемы в работе [8] была предложена конечно-объемная схема для решения двумерных уравнений Максвелла на треугольных сетках. Предложенная схема имела второй порядок аппроксимации и позволяла сохранить его даже для случая разрывной диэлектрической проницаемости. Однако построенные алгоритмы позволяли проводить расчеты только для случая, когда разрыв диэлектрической проницаемости проходил по координатной линии.

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

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

2. УРАВНЕНИЯ МАКСВЕЛЛА

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

ЭО Э?

■ - 10 Н = 0, В=еЕ,

ЭВ

-------+ 10 Е = 0.

Э?

= 0,

в = тн,

&уВ = 0,

где Е - электрическое поле, Н - магнитное поле, В - электрическая индукция, В - магнитная индукция, е - диэлектрическая проницаемость, ц -магнитная проницаемость. Далее везде полагаем, что магнитная проницаемость ц = 1. Данная система может быть представлена в векторной консервативной форме

- и + -Э?

Э -К + —К + — К=0,

Эx^

Эх,

(1)

Эх1 их2 их3

где и - вектор консервативных переменных, К1 К2 и К3 - векторы потоков.

и=

Г в ' Г 0 ^ Г - Н 3 ' Г Н 2 1

А Н 3 0 - Н1

Д В1 - Н 2 0 , К2 = Н1 Е3 , Е3= 0 - Е2

В2 - Е3 0 Е1

V В3 J V Е2 J V - Е1J V 0 J

Проинтегрировав уравнение (1) по объему П с границей ЭО, можно получить эквивален-тую интегральную форму

Э р

— Гшп +

Э?-!

| (п1К1 + п2 К2 + П3 К3 )

=0, (2)

где п = (п1, п2, п3) - внешняя нормаль. Систему (1) также можно записать в недивергентной форме

—V + Д—V + А2—V + А3—У = 0, (3)

Э? Эх1 Эх2 Эх3 ( )

где V - вектор потоковых переменных, связанных с вектором консервативных переменных

V =

0 =

Г Е 1

Е2

Е3

Н1

Н 2

VН 3 J

а матрицы А1, А2, А3 записываются как

V= ©и

Г1

0 0 0 0 0

е

0 1 0 0 0 0

е

0 0 1 0 0 0

е

0 0 0 1 0 0

0 0 0 0 1 0

0 0 0 0 0 1

Л

А1 =

А2 =

Г 0 0 0 0 0 0 ^ 1 е

0 0 0 0 0 1 е

0 0 0 0 0

0 0 0 0 0 0

0 0 -1 0 0 0

V0 Г 1 0 0 0 0 J 11

0 0 0 0 0 е

0 0 0 0 1 е 0 0

0 0 0 0 0

0 0 1 0 0 0

0 0 0 0 0 0

V -1 0 0 0 0 0 J

(

Aз =

0 0 0 0

1

0 0 0 — 0 0

000 0 -10

1 0 0

000

е

0

0

0

0

00

00

00

00

3. РАЗНОСТНАЯ СХЕМА

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

и барицентр XД как

оД = ш, хВ = — \хт.

А ОДд

Для каждой грани Г определим площадь 5Г

— г.

и центр Xе как

£г = , хС = — Гх^.

■, ^г *

г г г

Для приближенного решения уравнения (2) рассмотрим разностную схему

и п+1 и п т

X Г (п1Е1 + п2р2 + пзРз )^Г=0, (4)

-+^(п

к=1 г

где ип - аппроксимация значения и в барицентре 1-й ячейки хв в момент времени ?п = пт, т -шаг по времени, а ОД - объем / -й ячейки.

Для того чтобы в (4) найти значение и на новом временном слое - ип+', надо вычислить интегралы по граням Гк ячейки А/ которые представляют собой потоки искомых величин через грани. Предполагаем, что в ячейке искомые функции изменяются линейно, по значениям этих функций в барицентре ячейки и по вычисленным градиентам функций находим значения функций в центре грани со стороны /-й ячейки в момент времени пт + т/2. Аналогичным образом находим значения функций в центре грани со стороны ячейки, находящейся по другую сторону грани. По значениям функций по разные стороны грани, вообще говоря различным, находим потоки в центре грани в момент времени пт + т/2. Тогда интегралы в (4) приближенно вычисляем по формуле прямоугольников и получаем

5Д - площадь к-й грани К/ - поток через к-ю

/

грань.

3.1. Нахождение потоков через границу ячейки

Пусть п = (пь п2, п3) - вектор нормали к общей грани ячеек А^ и АЯ в точке Xе. Если предположить, что производные компонент электромагнитных полей по касательной к грани равны 0, то систему в недивергентной форме (3) можно переписать как

—V + А—V = 0,

Э? Эп

где А = А\п\ + А2п2 + А3п3. Для этой системы можно рассмотреть одномерную задачу Римана, где в качестве начальных значений принять аппроксимацию компонент электромагнитных полей в центре грани со стороны ячеек А^ и АЯ -У1(Хе) и УЯ(Хе). Решение этой задачи [7, 8] будем использовать для вычисления потока через границу ячейки.

Обозначим В = Я1АЯ, здесь В - диагональная матрица из собственных значений матрицы А, В± - диагональные матрицы, полученные из В заменой всех отрицательных (положительных) собственных чисел нулями, Я - матрица, столбцы которой являются правыми собственными векторами матрицы А. Тогда

А = ЯВЯ-1 = ЯВ+Я-1 + ЯВ-Я-1 = А+ + А~ .

Учитывая V = ©и, для вычисления потоков через грань в (5) будем использовать

К=©-1 (А+^ (Xе)+ А^Я (Xе))

или

Е=С+Уь (хС)+ С~УК (хС).

'1\х )~ге гЯ\х ). (6)

Матрицы е + и е задаются по формулам

Г 1 ^

—;=(Е - п п) л/е

N

-Ы е(п‘п - Е)

1

~^=(п*п - Е) - N

л/е

N е(Е - п‘п)

<

где

е = 2еьеК/(еь +ек),

' 0 - пз п2 "

N = пз 0 - п1

1- п2 п1 0 1

а Е - единичная матрица 3 на 3.

ип+\ = ип

а

(5)

Д. к=1

Т

3.2. Нахождение значений компонент электромагнитных полей на границе ячейки

Рассмотрим исходную систему уравнений в недивергентной форме (3). Пусть X - барицентр ячейки, а Xе - центр грани. Тогда ¥(Х) может быть найдена по следующей формуле со вторым порядком по времени и пространству:

У(Х е) = У(Х в) + —(Xв )(Х е

Эх

■X в )■ ЭУ

Л ^(Хв) + л2 ^(Хв) + Аз ЭУ(Хв)) Эх1 Эх2 Эх-

Л

Итак, значения электромагнитных полей на грани ячейки сетки в (6) выражаются следующим образом:

V (Xс)=У(Х В)+|^(Х В )(хс - X В) -

ох

/г ( ОУ ^ Злл Л

ЭУ

ЭУ

Лі - (х в) + Л2—(Х в) + Лз—(Х в)

Эх1

Эх

Є\ — лг/^вх , ЭУ /V^уе

Уд (Xе) = у(х в) +—(X д )(Х

Эх

Лі |^(Х в) + Л2 |^(Х д) + Лз^(Х д)

Эх3

X д )■

ЭУ

л

Эх1

Эх

Эх

3.3. Нахождение градиентов компонент электромагнитных полей в ячейке

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

( Е Л

Ж = 5(0, ф, е)

Е2

V Ез у

^ 88Іп(0)сО8(фф 88Іп(0)8Іп(фф 8008(0) ^ СО8(0О8(0)ф) СО8(0О8(0)ф) - БІп(0І

- 8Іп(фІ СО8(фО

0

где 0 и ф - углы вектора нормали к поверхности разрыва в сферической системе координат, а 8 -диэлектрическая проницаемость в ячейке. В остальных ячейках выберем Ж = Е. В узлах сетки

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

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

|Х ;, Е(0,, фг, е}.)Е}.^

где j пробегает индексы ячеек, которые используются для вычисления градиентов в ячейке /. После нахождения градиентов в /-й ячейке сетки по найденным градиентам вычислим значения Ж/ (Хр )=Ж/Р в ее узлах Хр={хр }

Ж/ (хр) = 3(0р, фр, 8/ )Н-1(0/, ф/, 8/ )

Ж + -°-Ж (ХВ )(хр - хВ

V Ох1

+ О^Ж/ (ХВ )(хр - х2В/ )н

+

+

Э в

+ ЭкЩ(X ; )(х

(хзР - хзвг)

где 0р и фр - значения углов в узле р. В одном узле X" получаем несколько различных значений Ж/(Хр) по числу соседних ячеек, за итоговое значение Ж = Ж(Х) принимаем их среднее арифметическое. Теперь в каждой ячейке найдем градиенты Е с помощью метода наименьших квадратов. Для этого будем использовать значения Е в вершинах ячейки. Таким образом, в методе наименьших квадратов для вычисления градиентов Е в ячейке / возьмем

Xрк, 5-1(0рк

кк Р ,е,■ )Жр )

где рк, к = 1,2,... - номера вершин ячейки /.

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

5. ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ

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

декомпозиции сеток позволяет разбивать на фрагменты сетки, построенные с помощью Ош8Ь [10], ТЕТСЕК [11] и КЕТСЕК [12]. Каждый процесс вычислительной программы проводит расчеты на своем фрагменте. Для передачи данных между процессами используется интерфейс передачи сообщений МР1.

5.1. Декомпозиция сетки

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

5.2. Вычислительная программа

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

• передача значений электромагнитных полей в ячейках;

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

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

• вычисление окончательных градиентов;

• вычисление значений на гранях;

• передача значений на гранях;

• вычисление значений в ячейках на новом временном слое.

Особо отметим, что алгоритм не требует обмена градиентами между процессами.

6. РЕЗУЛЬТАТЫ ТЕСТОВЫХ РАСЧЕТОВ

Для проверки свойств предложенной схемы были проведены тестовые расчеты, в которых использовались тетраэдральные сетки. Для построения тетраэдральных сеток использовался открытый программный продукт Ош8Ь [10]. Точность численного алгоритма оценивалась путем сравнения с аналитическими решениями. Ошибка численного решения в момент времени (п = пт вычислялась по формуле

\\у" (хв)-Vехас‘ (Хв,гп )

2 _

V ехас‘ (хв, Ґ)

Т £ і=1 £ V (хв')-ГГ' ХВі, 1п)) к=1 • 5 д. і

Т 1? £ (г™ (хВі, 1п)) _ к=1 _ • 5 А

где Т - общее количество тетраэдров в вычис-

лительной области, Г^уХ ’ I и Гкехас \Х ' ,Ґ ) -вычисленные и точные значения электромагнитных полей в барицентре ячейки і, соответственно.

6.1. Тест 1

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

г = уІ х1 + х22 = а,

е = е(г) = <

£1 = п , 0 < г < а,

е2 = п2, г > а.

В этом случае система уравнений Максвелла имеет аналитическое решение [13], г - компоненты элекромагнитных полей которого в цилиндрических координатах записываются в виде

Ег = У1 ^ иг | • ео8(0)со8(к - —г),

Нг = -—&1л иг | • 8Іп(0)ео8(кґ - —г),

к І а I

ь

Нг = --— sJ11 иг | • 8Іп(0)со8(кґ - —г),

к ^ а )

при г < а, и

Е = -(и)

к1 -—г I • со8(0)со8(к? - —г),

КД—) {а )

Н, = --— 5 -1(и) К1 (—г1 • 8Іп(0)со8(к? - —г),

к К1 ^ а

при г > а, где ^ и - функции Бесселя первого рода, К0 и К1 - функции Бесселя второго рода,

50 = 5—2/к2П02 ,

2 2

5 = 5—2/к2п2, и = а^к2п12 -—2 , л/к2П22 -—2 ,

“|-1

-о(и) - - 2 (и) К о(—) - К 2 (—)

,,22 п2

— = ^/к п -

и-1 (и)

—К1 (—)

а в находится из дисперсионного соотношения

-1(и). + . К1'(—)

и-1 (и) —К1 (—)

+

1 1

"Г + ~ и —

-1(и) , по2 К1'(—) и-1 (и)

по2 1 "

—К1 (—)

2 + 2 2 п—

Остальные компоненты электромагнитных полей могут быть выражены через г-компонен-ты с помощью уравнений Максвелла.

Тестовые константы 81 = 2,25, 81 = 1,0, к = = 6,0, а = 0,64, в = 8,402440923258, и = = 2,063837416842, = 3,764648073438.

В качестве расчетной области был выбран цилиндр радиуса 1.0 и высотой 0.6:

^х2 + х2 < 1,0, 0 < х3 < 0,6.

Расчеты проводились на последовательности сеток, состоящих из 18845, 51147, 144048, 405279 и 1250790 тетраэдров. Сетки строились таким образом, чтобы разрыв диэлектрической проницаемости проходил по граням тетраэдров. Шаг по времени выбирался пропорциональным линейным размерам тетраэдров. На рис. 1 показан пример расчетной сетки, состоящей из 18845 тетраэдров. Не показаны тетраэдры, у которых первые две координаты барицентра меньше 0,5. Подобласть с более высокой диэлектрической проницаемостью 81 = 2,25 обозначена темным цветом. На рис. 2 показаны изолинии распределения третьей компоненты

магнитного поля Н3 в момент времени Т = 4,73 в сечении х3 = 0,3, полученные на сетке из 144048 тетраэдров. На рис. 3 показана эволюция ошибки 52 в норме Ь2 на последовательности из пяти сеток. В табл. 1 приводятся максимальные значения ошибки в норме Ь2. Как следует из результатов расчетов, алгоритм позволяет проводить расчеты со вторым порядком точности в областях с разрывной диэлектрической проницаемостью в случае, когда разрыв проходит по криволинейной поверхности.

Рис. 1. Пример расчетной сетки

Рис. 2. Магнитное поле г-компонента

Таблица 1

Максимальная ошибка

Число тетраэдров §2 Порядок аппроксимации

18845 0,055064

51147 0,027107 2,12

144048 0,013602 2,06

405279 0,006120 2,14

1250790 0,002880 2,11

Рис. 3. Эволюция ошибки

Рис. 4. Электрическое поле х компонента 6.1. Тест 2

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

e = e(r) =

I e = n2, 0 < r < a,

I e2 = n2, r > a.

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

Расчеты проводились с использованием последовательной и параллельной реализации предложенного алгоритма на кластере Новосибирского государственного университета. Кластер построен на базе блэйд-серверов HP BL460c, имеющих по 16 ГБ оперативной памяти и по два 4-ядерных процессора Xeon 5355 (2,66 ГГц). В качестве коммуникационной среды использовался InfiniBand. В параллельной реализации использовалась декомпозиция расчетной области на 2, 4, 8, 16, 32, 64, 128 частей. Результаты всех расчетов с использованием параллельной версии и расчета последовательной версии совпадали. На рис. 4 показаны изолинии распределения первой компоненты электрического поля E1 в момент времени T = 2,00 в сечении х1 = 0,0. Время счета было различным. В табл. 2 приводятся затраты времени на проведение расчета в зависимости от числа процессов. Видно, что достигается ускорение близкое к линейному, что говорит о высокой эффективности параллельной реализации.

Т аблица 2

Ускорение

Число процессов Время счета, с

1 90179

2 48512

4 25100

8 15243

16 7052

32 3333

64 1715

128 939

ВЫВОДЫ

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

5

2

T

Волна распространяется в направлении х3 и имеет компоненты электромагнитных полей Е\ и И2. Для этой задачи система уравнений Максвелла имеет аналитическое решение [14], но в силу его громоздкости приводить его здесь не будем.

СПИСОК ЛИТЕРАТУРЫ

1. Yee K. S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // IEEE Trans. Antennas Propagat. 1966. V. 14. P. 585-589.

2. Taflove A. Advances in Computational Electrodynamics: the Finite-Difference Time-Domain Method // Boston, Artech House, 1998.

3. Taflove A. and Hagness S. C. Computational Electrodynamics: the Finite-Difference Time-Domain Method. Boston, Artech House, 2000.

4. Sullivan D. M. Electromagnetic Simulation Using the Finite-Difference Time-Domain Method // New York, IEEE, 2000.

5. Hermeline F. Two coupled particle-finite volume methods using Dalaunay-Voronoi meshes for ap-plorimation of Vlasov - Poisson and Vlasov - Maxwell equations // J. Comput. Phys. 1993. V. 106. P. 1-18.

6. Cioni J.-P., Fzoui L., Issautier D. Higher order upwind schemes for solving time domain Maxwell equations // La Recherche Aerospatiale. 1994. №. 5 P. 319328.

7. Лебедев А. С., Федорук М. П., Штыри-

на О. В. Решение нестационарных уравнений Максвелла для сред с неоднородными свойствами методом конечных объемов // Вычисл. технологии. 2005. Т. 10, № 2. С. 60-73.

8. Лебедев А. С., Федорук М. П., Штыри-на О. В. Конечно-объемный алгоритм решения нестационарных уравнений Максвелла на неструктурированной сетке // ЖВМ и МФ. 2006. Т. 47, №. 7. C. 1286-1301.

9. Cioni J.-P., Fzoui L., Steve H. A parallel time-domain Maxwell solver using upwind schemes and tri-

angular meshes // IMPACT Comput. Sci. Eng Academic Press, Orlando, FL, USA 1994. V. 5. P. 215-247.

10. Geuzaine C. and Remacle J.-F. Gmsh: a threedimensional finite element mesh generator with built-in pre- and post-processing facilities// Int. J. Numer. Meth. Engng. 2009. V. 3 P. 1-24.

11. Hang Si. Adaptive tetrahedral mesh generation by constrained Delaunay refinement // Internat. J. Nu-mer. Methods Engrg. 2008. V. 75. P. 856-880.

12. Schoberl J. NETGEN - An advancing front 2D/3D-mesh generator based on abstract rules // Com-put. Visual. Sci. 1997. V. 1. P. 41-52.

13. Okamoto K. Fundamentals of Optical Waveguides. London, Academic Press, 2000.

14. Bohren C. F., Huffman D. R. Absorption and scattering of light by small particles, Wiley, 1998.

ОБ АВТОРЕ

Исмагилов Тимур Зинферо-вич, ст. преп. Дипл. бакалавр по прикл. математике (НГУ, 1998). Дипл. магистр по математике (НГУ, 2000). Иссл. в обл. вычислительной математики, математического моделирования, параллельных численных методов.