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

УДК 519.642

Д. А. Миронов

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

Аннотация. Рассматривается задача дифракции стороннего электромагнитного поля на локально неоднородном теле, помещенном в свободном пространстве. Поставленная задача сводится к объемному сингулярному интегральному уравнению. Решение задачи производится параллельно численным методом Галеркина и численным методом коллокации. В связи с большой емкостью решение задачи численным методом Галеркина при различных параметрах было реализовано с использованием двух программных продуктов для супер-компьютерных вычислительных комплексов: реализации MPI и программной системы x-com. Исследованы особенности выполнения задачи на суперкомпь-ютерном комплексе.

Ключевые слова: дифракция стороннего электромагнитного поля, объемное сингулярное интегральное уравнение, метод Галеркина, метод коллокации.

Abstract. In this paper the problem of diffraction of external electromagnetic field on locally heterogeneous body placed in free space is considered. The formulated problem is reduced to three-dimensional singular integral equation. The problem is solved using Galerkin numerical method and at the same time using the collocation numerical method. Because of high capacity the task solving was realized by Galerkin numerical method at various parameters with use of two types of software for supercomputing complexes: realization of MPI and realization of program system x-com. The features of performance of a task on a supercomputer complex are investigated.

Keywords: diffraction of external electromagnetic field, three-dimensional singular integral equation, Galerkin method, collocation method.

Постановка задачи для системы уравнений Максвелла

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

магнитной проницаемостью ^q и положительной (3 х 3)-матрицей-функцией (тензором) диэлектрической проницаемости ё(х). Компоненты ё(х) являются ограниченными функциями в области Q , e е Lx (Q), а также ё 1 е Lx (Q). Граница dQ области Q кусочно-гладкая.

Требуется определить электромагнитное поле E, H е Lq (Q), возбуждаемое сторонним полем с временной зависимостью вида e-l(at. Источник стороннего поля - электрический ток j0 или падающая плоская волна.

Будем искать электромагнитное поле E, H, удовлетворяющее уравнениям Максвелла, условиям непрерывности касательных компонент поля при переходе через границу тела и условиям излучения на бесконечности:

rotH = -ШёЁ + jE ; rotE = Iro^QH; (1)

1E 1| ае =[ H ]1 ае = °;

ае

( E Л ( E Л

дЯ

v H J

— ik°

v H J

( E Л

= o(R—1), = O(R—!), R := Ixl

v H J

(2)

(3)

Здесь ко - волновое число свободного пространства (вне Q ).

Краевую задачу (1)-(3) можно свести к объемному (векторному) сингулярному интегральному уравнению [1]:

E (х) +

g ( х)

g°

—I

E (x) — v.p.j Гі( x, y)

е

g( y)

g°

—I

E(y)dy —

j Г( x, y)

е

g( y)

g°

—I

E(y)dy = E (x),

(4)

где

Г( x, y) = k°G(r) + (•, grad) grad G°(r); Гі( x, y) = (•, grad) grad Gi(r). Функция Грина имеет вид

1 gik°lx—yl

G( x, y) =

4л | x — y |

G(r) = G°(r) + Gi(r), r =| x — y |; G°(r) =

eikor — і 4лг

Gi(r) =

4nr

Численные методы решения

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

По методу Галеркина решение интегрального уравнения сводится к решению системы линейных алгебраических уравнений (СЛАУ) вида [1]:

AX = B,

(5)

где

(A11 A12 A13 Л ( B1 Л

A = A21 A22 A23 , b = B2

V A31 A32 A33 J V B3 J

Ак1 - блок матрицы вида:

Akl = j ^klfj (x )fik (x )x — sklk° jj G (x, y )fj (y )fik (x )dydx +

nj nnk

nk nj

а

П кг П1].

+1 I сУ(у)^Т/к (х^х;

д_

дхк

Б1к = (Ерк,/к), к,I = 1,2,3; /, ] = 1,...,

/к1т

1 - “11 х1 - х1,к (, хе П1/т,

О, х ^ П1/т ;

П11т {х : х1,к-1 ^ х1 ^ х1,к+1, х2,1 ^ х1 ^ х2,1+1, х3,т ^ х3 ^ х3,т+1};

а2 - а1 , Ь2 - Ь , . с2 - с

х1к = а1 + —------------ к, х21 = Ь1 + 2—--------- I, хт,т = С1 + 2—--------- т,

п п п

к = 1,..., п -1; I, т = 1,..., п/2 -1;

(6)

к1 :=| х1,к - х1,к-11.

2 3

Функции /к1т , /к1т, зависящие от переменных х2 и х3 соответственно, определяются аналогичными соотношениями.

По методу коллокации решение интегрального уравнения сводится к решению системы линейных алгебраических уравнений (СЛАУ) вида [2]:

3 1

(А» р1 ■ I Р р1т^рт -г-II Брд1т^дт Ер1,1 1,..., 3,

т=1 £р q т=1

(7)

где

Б

pqln

3 2к0 ,2 ^(хр1 - х1)(хрп - хп)

2 - кО

г2 г

+1 кр+*р - -1

§1п

Б

'рр1п ®1 §1п + §1п | | к0 С(г)

ёх^2йх^, р Ф q; 1 -(хр^х^хрп2х2

кр2

4лг

3Ф(г)(хр1 - х1)(хрп - хп )

-Ф (г)

>йх1<Лх2йх3, р = q;

г =

1 (хп хрп )" ,п=1

1/2

а1 = / (к\, к2, кэ), а2 = / (/¡2, кь Л3), а3 = / (¿3, кь /¡2);

f Й Ä2, Йз) =1----------*

к

arcsin

йЙз

VhF + hi

+ arcsin

yjh2 + h2

г| + h32

Ф(г ) = 1 + (1 - ik§r )

всюду дифференцируемая функция;

exp(ik0r) - 1 - ikor (kor)2

,fa И2 Из

xpi = M + Y xp2 = P2И2 + ^ xp3 = Р3И3 + Y

î Й И2 И3

xq\ = qihL + Y xq2 = Ч2И2 + ^ xq3 = Ч3И3 + y =

Xp/ -1 -я декартова координата p -й узловой точки;

тт Щ Щ п2 п2

П p = \ x : xp1 - Y < x1 < xp1 + Y xp 2 2 < x2 < Хр 2 2 ,

m m

xp3 - — < Х3 < xp3 + —

т-г I Й Й Й h2

^ = \ x : xqi - у < x1 < xqi + ^, xq2 -у < x2 < xq2 + ^,

"3 ^3 I

Xq3 < X3 < Xq3 + y

p = (pi,p2,p3), pi = °,--,N1 ~ 1, p2 = °,--,N2 -1, p3 = 0,...,N3 -1; q = (qi,q2,q3), qi = °,-..,Ni -1 q2 = °,---,N2 -1, q3 = °,...,N3 -1.

Уравнения (5) и (7) решаются методом сопряженных градиентов [3] -итеративный метод вида

X* = A • X*-1,

где X* - решение уравнения на i -й итерации; X° = B .

Итерации выполняются до тех пор, пока не будет выполнено условие

| X1 - Xi-i |< е, где е - заданная точность.

Для вычисления коэффициентов Aj были задействованы ресурсы многопроцессорных вычислительных комплексов с использованием двух программных интерфейсов - MPI ( m < 8) и x-com ( m = 16). Для вычисления Bpq/n использован программный интерфейс MPI. Далее будут освещены особенности реализации алгоритмов с использованием данных программных интерфейсов.

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

Для уменьшения времени на работу алгоритма программы и уменьшения объема занимаемой памяти в [4] предложено учитывать симметрию коэффициентов АН (6) - достаточно вычислить и хранить в памяти коэффициенты блоков Ац и Ау1 . Общее количество коэффициентов данных блоков

вычислялось по формуле

N = (т + 1)6 • 2, (9)

где т - количество интервалов разбиения всей области по одной координате.

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

N = т6 • 2. (10)

Зависимость количества коэффициентов АЦ и необходимого объема

оперативной памяти для их хранения прямо пропорциональная. Результаты расчета необходимого объема оперативной памяти при различных значениях т по формулам (9) и (10) представлены в табл. 1.

Таблица 1

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

т Количество Мбайт при использовании формулы (9) Количество Мбайт при использовании формулы (10)

4 0,477 0,125

8 16,218 8

16 736,620 512

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

Таблица 2

Коэффициенты уменьшения использованного объема оперативной памяти

т Коэффициент уменьшения при использовании формулы (9) Коэффициент уменьшения при использовании формулы (10)

4 4,5 17,17

8 4,5 9,123

16 4,5 6,474

Количество коэффициентов Врц1п (8) при N1 = N2 = N3 = т вычисляется формуле

N = т6 • 9. (11)

В работе [2] показано, что для решения интегрального уравнения достаточно вычислить малую часть коэффициентов Врц1п, от которых зависят

все коэффициенты для уравнения (7) (последние можно быстро вычислить в дальнейшем). Количество хранимых коэффициентов Bpqln вычисляется по

формуле

N=т3•6. (12)

В табл. 3 приведены необходимые объемы памяти для хранения коэффициентов матрицы уравнения (7) и необходимых коэффициентов Bpqln .

Таблица 3

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

m Количество Мбайт при использовании формулы (11) Количество Мбайт при использовании формулы (12)

4 0,5625 0,00586

8 36 0,04688

16 2304 0,375

Реализация MPI-версии алгоритма вычисления коэффициентов. Алгоритм распределения вычислений на многопроцессорных комплексах

MPI [5] - удобный стандартный API для использования в прикладных задачах ресурсов многопроцессорных комплексов. На каждом вычислительном многопроцессорном комплексе используется одна или несколько реализаций (компиляторов) MPI.

Для упрощения вычислений необходимых коэффициентов Aj (6), Bpqin (8) и передачи их между процессорами на многопроцессорных комплексах память выделяется в виде одномерного массива.

Общая схема алгоритма вычисления коэффициентов с учетом использования MPI в реализации [3] представлена на рис. 1. Количество выделенных процессоров на задачу - p .

Количество коэффициентов C , которое необходимо вычислить на каждом процессоре, вычисляется по формуле

C =

N

p.

N

p

І NІ

+1, если номер процессора меньше < — !>,

ІР \

, если номер процесора больше или равен

N

где N - общее количество коэффициентов; { } - остаток от целочисленного деления; [ ] - целая часть деления; р - количество выделенных процессоров на задачу.

1. Вычисление коэффициентов

Массив

1 г 1 г коэффициентов і г

Вычисления Вычисления Вычисления

на процессоре 0 на процессоре 1 на процессоре

(Р - 1)

запись результатов на процессоре 0 и выход

Рис. 1 Общая схема алгоритма вычисления коэффициентов с использованием многопроцессорных комплексов

Программы, реализованные с использованием МР1-функций, были запушены на суперкомпьютерном комплексе СКИФ МГУ. Основные характеристики комплекса представлены в табл. 4.

Таблица 4

Основные характеристики суперкомпьютерного вычислительного комплекса СКИФ МГУ

Модель процессора Количество процессоров Минимальный объем оперативной памяти на один процессор

Intel Xeon E5472 3.0 ГГц от 1 до 5000 в зависимости от количества и объема задач, уже работающих на комплексе 2 Гбайт

Для вычисления коэффициентов блоков Лц и Л-12 , при т = 8, п = 10, произведен запуск программы. Количество выделенных процессоров на задачу - 500. В табл. 5 показано время выполнения программы. Для сравнения приведено время выполнения представленного в [4] алгоритма и алгоритма с учетом выявленного свойства, т.е. при количестве коэффициентов, вычисленных по формулам (9) и (10) соответственно.

Таблица 5

Время на вычисление коэффициентов в блоках Лц и Л12

Время на вычисление элементов матрицы При количестве коэффициентов, вычисленном по формуле (9) При количестве коэффициентов, вычисленном по формуле (10)

Секунды 6012,26 2241,74

Минуты 100,204333333 37,362333333

Часы 1,670072222 0,622705555

По результатам, представленным в табл. 5, нетрудно вычислить примерное время выполнения вычисления коэффициентов матрицы при т = 16, п = 10. При п = 500 необходимо 40-45 ч на выполнение с использованием алгоритма, учитывающего выявленное свойство.

Вычисление коэффициентов матрицы по методу Галеркина.

Алгоритм распределения вычислений при помощи системы «х-сот»

Основные причины применения системы «х-сот» для решения задач.

Основные понятия системы «х-сот»

При заполнении блоков Лц и Л12 на суперкомпьютерам вычислительном комплексе СКИФ МГУ с использованием МР1-функций при т = 16, п = 10 с учетом дополнительных расходов системы на передачу коэффициентов между процессорами, время на решение задачи составит ~ 1-2 суток при 1000 выделенных процессах. Если в течение этого времени возникнет ситуация, при которой хотя бы один процесс прекратит работу, необходимые коэффициенты блоков мы не получим.

Программная система «х-сош», разработанная специалистами НИВЦ МГУ, предназначена для многопроцессорных комплексов. Основная цель системы - решать задачи в многопроцессорных средах, где возможно прекращение работы одного или нескольких процессоров во время решения задачи. Эта система оптимальна для решения задач, которые допускают разбиение на независимые подзадачи.

Основные термины системы «х-сош»:

1. Порция - независимая подзадача общей задачи. Может выполняться одновременно, параллельно с другими подзадачами (порциями).

2. Серверная часть системы «х-сош». Программный модуль, содержащий алгоритмы распределения порций по процессорам. Содержит функции, определяющие:

а) алгоритмы разбиения всей задачи на порции - количество порций, время на выполнение одной порции. Эти функции специфичны для конкретной задачи;

б) алгоритмы объединения результатов работы всех подзадач-порций.

3. Клиентская часть системы «х-сош». Программный модуль, содержащий алгоритмы приема на выполнение порции каждым процессором. Выполняет алгоритм работы подзадачи.

Для нашей задачи при работе в системе «х-сош» необходимо реализовать:

1. Функции серверной части, реализующие алгоритм разбиения всей задачи на порции.

2. Алгоритм работы каждой порции.

Функции серверной части, реализующие алгоритм разбиения всей задачи на порции

После выполнения МР1-программы алгоритма вычисления коэффициентов блоков Лц и Л12 вычислено среднее время на вычисление одного коэффициента - ¿0.

Задано среднее время выполнения порции - ¿1.

Известно N - общее количество коэффициентов блоков Лц и Л12 • Время T для вычисления коэффициентов блоков Лц и Л12 определяется по формуле

T = to • N •

Тогда количество порций P определяется по формуле

T

P =

где [ ] - целая часть деления.

Количество вычисляемых коэффициентов в подзадаче (порции)

В каждой порции выполняется вычисление определенного количества коэффициентов.

Перед выполнением алгоритма решения подзадачи, клиенту передаются:

- номер порции;

- количество порций.

Количество коэффициентов, вычисляемое в порции с номером 7 = 1, Р , определяется по формуле

С/ =

N

P

N

Т

+ 1, если номер порции меньше или равно | р если номер порции больше | Р

где N - общее количество коэффициентов в блоках Лц и Л12 , Р - количество порций, { } - остаток от целочисленного деления, [ ] - целая часть деления.

Статистика по работе системы «х-сот» при выполнении задачи вычисления коэффициентов

Программа с реализацией алгоритма вычисления коэффициентов блоков Лц и Л12 при т = 16, п = 10 была запущена в системе «х-сош». Было задействовано 4000 процессоров.

Время счета на «х-сош»: общее время счета: 06 ч 54 мин 23 с (24863 с).

Общее количество коэффициентов в блоках Лц иЛ12 : 33554432.

Объем памяти для хранения коэффициентов блоков Лц и Л12 : 536870912 Байт (524288 кБ; 512 МБ; 0,5 ГБ).

Время для расчета с использованием только одного процессора для запуска программы с использованием «х-сош»: 19525 ч 55 мин 33 с (813,54 сут; 2,23 лет).

Количество частей-порций в задаче, которые могли выполняться в произвольном порядке (одновременно, параллельно): 140428.

Количество коэффициентов в одной части-порции: от 199 до 200.

Время расчета каждой порции:

- минимальное: 458,522 с (7,64 мин);

- максимальное: 609 с (10,15 мин);

- среднее: 485,211 с (8,09 мин).

Количество отправленных порций на выполнение (порция может отправляться более одного раза): 168062. Из них посчитанных порций (порция может считаться более одного раза): 144875.

Количество отправок без посчитанных порций: 23187.

Количество отправок с лишним счетом порций: 4447.

Данные статистики выполнения алгоритма на системе «x-com»:

1. Доля эффективности системы только с учетом проблем на узлах (не-посчитанные порции) относительно идеальной ситуации, когда каждая порция будет посчитана с первого раза и только один раз, составляет 0,86203 (86,203 %);

2. Доля эффективности системы с учетом дополнительных расходов на качественное выполнение задачи (лишний счет порций) относительно идеальной ситуации, когда каждая порция будет посчитана с первого раза и только один раз, составляет 0,96930 (96,930 %).

3. Доля эффективности системы с учетом проблем на узлах и дополнительных расходов на качественное выполнение задачи относительно идеальной ситуации, когда каждая порция будет посчитана с первого раза и только один раз, составляет 0,83557 (83,557 %).

Важный результат при применении системы «x-com»

Применение системы «x-com» позволило получить коэффициенты блоков Ац и Ayi при m = 16, n = 10. Ранее это не удалось сделать по следующим причинам:

1) большая вычислительная сложность задачи;

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

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

Для вычисления необходимых коэффициентов Bpqin произведен запуск программы при m = 16, n = 10. Количество выделенных процессоров на задачу - от 1 до 3. В табл. 6 показано время выполнения программы.

Таблица 6

Время вычисления коэффициентов Bpqin

Один процессор Два процессора Три процессора

Время на вычисление элементов матрицы, с 8,033112 4,018449 2,677404

Время на вычисление элементов матрицы, включая запись на жесткий диск результатов, с 12 11 8

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

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

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

Действительно, по скорости вычисления коэффициентов метод Галер-кина явно уступает методу коллокаций. Но необходимо отметить, что при реализации алгоритма по методу Галеркина получены не только коэффициенты, но и результаты решения задачи в целом после решения СЛАУ при параметрах m = 16, n = 10 (впервые). Кроме того, испытана при вычислении коэффициентов система «x-com», показавшая высокое качество работы над трудоемкой с вычислительной точки зрения задачей. По сравнению с применением среды MPI, среда «x-com» предназначена для многопроцессорных вычислительных систем, где возможно прекращение работы нескольких процессоров во время вычислений.

Список литературы

1. Медведик, М. Ю. Применение ГРИД-технологий для решения объемного сингулярного интегрального уравнения для задачи дифракции на диэлектрическом теле субиерархическим методом / М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2008. - № 2 (6). - С. 2-14.

2. Самохин, A. Б. Интегральные уравнения и итерационные методы в электромагнитном рассеянии / A. Б. Самохин. - М. : Радио и Связь, 1998.

3. Ор тега, Дж. Введение в параллельные и векторные методы решения линейных систем / Дж. Ортега. - М. : Мир, 1991.

4. Миронов, Д. А. Применение суперкомпьютерных вычислительных комплексов для решения объемного сингулярного интегрального уравнения задачи дифракции на диэлектрическом теле / Д. А. Миронов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2008. - № 3. - С. 55-62.

5. MPI: A Message - Passing Interface Standart. Versión 1.0. - University of Tennessee, 1994. - May, 5.

Миронов Денис Алексеевич Mironov Denis Alexeevich

аспирант, Пензенский государственный Postgraduate student,

университет Penza State University

E-mail: dionis-m@mail.ru

УДК 519.642 Алехина, М. А.

Применение суперкомпьютерных вычислительных сред для решения объемного сингулярного интегрального уравнения задачи дифракции на диэлектрическом теле / Д. А. Миронов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2009. -№ 2 (10). - С. 14-24.