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

УДК 517.96

Д. А. Миронов

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

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

Введение

В работе рассматривается задача дифракции стороннего электромагнитного поля на локально неоднородном теле, помещенном в свободном пространстве. Актуальность работы определяется применением результатов исследования, например, при решении задач дифракции в СВЧ-диапазоне. Для численного решения задачи использован метод объемных сингулярных интегральных уравнений, актуальность использования которого обоснована в [1].

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

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

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

Требуется определить электромагнитное поле Е, Н е ^), возбуждаемое сторонним полем с временной зависимостью вида е-г(°1. Источник стороннего поля - электрический ток у0 или падающая плоская волна.

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

гоН = -1юё.Е + уЕ; ю\Ё = /юцоН; (1)

[ Ё]1 ю=[ Н С=°; (2)

_д_

дЯ

( Е Л ( Е Л -1 ( Е Л

Н

-

Н

= о(Я-1), = 0(Я-1), Я := |х| (3)

Н

V У V У V У

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

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

ё( x)

£o

-1

E (x) - v. p.J Гг( x, y)

Q

ё( У)

eQ

-1

E (y )dy -

J Г( x, y)

Q

£( У)

£0

-1

E(y)dy = E (x),

(4)

где

Г( x, y) = k0G(r) + (• ,grad)grad Go(r), Гг( x, y) = (•, grad) grad Gl(r). Функция Грина имеет вид

1 eiko|x-У\

G (x, y) = 4 \ \

4я | x - y |

G(r) = Go(r) + Gl(r),r =| x - y I; Go(r) =

eikQr -1 4nr

Gl(r)=

4nr

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

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

AX = B,

(5)

где

' All A12 A13 '

A = A21 A22 A23 , B = B2

1A31 A32 A33 J 1B3 J

где Акі - блок-матрицы вида

Аы = | (х)Ак (х)*- 8^0 | | G(хуУ] (у)$ (x

П j пП

J 1

+ J J G(x,у)д—fjl (у)д—fi (x)dydx;

Пk П J dyl dxk

1 J

Bk = (Eq , fk), k, l = 1,2,3; i, J = 1, ..., N,

fklm ='

(б)

1 - тг1 xl- xl,kl, x Є Пкт,

Q, x І Пkrlm ;

Пк1т {х : х1,к-1 < Х1 < х1,к+1, х2,1 < Х1 < х2,1+1, х3,т < х3 < х3,т+1};

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

Х1к = а1 +----- к, Х21 = «1 + 2—----1, Х3 к = С1 + 2—--- ш,

и и и

где к = 1,..., и -1; I, т = 1,..., и/2 -1; й1 :=| Х1 к - х^ к—1 I; и - количество

интервалов интегрирования по каждой координате.

2 3

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

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

X1 = А • Xм,

где X1 - решение уравнения на /-й итерации,

X0 = В.

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

IX' - Xм |<е,

где £ - заданная точность.

Учет симметрии. Распределение вычислений.

Алгоритм работы программы на многопроцессорных комплексах

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

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

N = (ш +1)6 2,

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

Для упрощения передачи данных между процессорами все коэффициенты блоков Ац и А12 матрицы хранятся в одномерном массиве {АI}, где индекс

I = (т + 1)( т +1)( т +1)( т +1)( т + 1)2/1 + (т +1)( т +1)( т + 1)( т +1)2/ 2 +

+(т + 1)(т + 1)(т +1)2/3 + (т + 1)(т +1)2 ]1 + (т +1)2 ] 2 + 2 ]3 +1,

где /1,/2,/3 - 0...т - составляющие индекса / по каждой из трех координат; Д]2,]3 - 0...т - составляющие индекса ] по каждой из трех координат

Допустим, нам доступно процессоров р > 1. Нумерация процессоров с нуля. Количество коэффициентов С вектора {А1} , которое необходимо вычислить на каждом процессоре, стараемся распределить равномерно:

С =

N

р.

N

р

I NI

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

IР)

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

где N - общее количество коэффициентов блоков Ац и Ац матрицы; {И/р} -

>"

остаток от целочисленного деления;

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

Общий объем векторов, необходимых для решения СЛАУ (5) (не включая вектор {} ), вычисляется по формуле

N0 = (т + 1)33 • 4.

Для решения СЛАУ (5) все значения элементов векторов, включая элементы вектора {ЛI} , необходимы на каждом процессе.

Общий объем векторов, включая вектор {ЛI}, необходимо учитывать, т.к. объем оперативной памяти любого ПЭВМ или суперкомпьютерного вычислительного комплекса всегда ограничен. Объем занимаемого места элементами векторов зависит от значения т .

Пользователю предоставляется возможность выбора значения т из множества {4,8,16,32} - выбранное значение зависит от субъективного выбора пользователем желаемого количества значений вектора решений в области.

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

Таблица 1

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

т Количество элементов матрицы с учетом симметрии Количество байт Количество кбайт Количество Мбайт Количество Гбайт

4 31250 500000 488,2813 0,476837 0,000466

8 1062882 17006112 16607,53 16,21829 0,015838

16 48275138 7,72х108 754299 736,6201 0,719356

32 2,58х109 4,13х1010 40358374 39412,47 38,48874

Таблица 2

Результаты расчета необходимого объема для хранения элементов векторов, необходимых для решения СЛАУ (без учета элементов матрицы)

т Общий размер Количество байт Количество кбайт Количество Мбайт Количество Гбайт

4 1500 24000 23,4375 0,022888 2,24 х10-5

8 8748 139968 136,6875 0,133484 0,00013

16 58956 943296 921,1875 0,899597 0,000879

32 431244 6899904 6738,188 6,580261 0,006426

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

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

Вектор матрицы

1 1 1

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

на процессоре 0 на процессоре 1 на процессоре (р-1)

2. Решение СЛАУ Одна итерация решения СЛАУ

Вектор решения

I Вычисления на процессоре 0 1 Вычисления на процессоре І 1 Вычисления на процесссоре (p-i)

Обмен данными

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

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

Программа была запушена на суперкомпьютерном комплексе СКИФ МГУ. Основные характеристики комплекса представлены в табл. 3 [3].

Таблица 3

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

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

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

Также на данном комплексе доступна возможность использования в программах МР1-функций для распределения вычислений и передачи данных между используемыми процессорами [4].

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

Основное требование при запуске на кластере комплекса СКИФ МГУ -определить до запуска время расчета, когда программа гарантированно завершит работу [6].

С учетом доступного объема оперативной памяти (см. табл. 3) и рассчитанного значения необходимого объема для хранения всех векторов (см. табл. 1, 2) выбранное значение т не должно превышать 16.

Для определения времени расчета при т = 8 произведены модельные запуски программы при т = 4, п = 3 с разными количествами используемых процессоров. Данные о времени выполнения программы при модельных запусках представлены в табл. 4

Таблица 4

Данные о времени выполнения программы при модельных запусках (т = 4, п = 3)

Количество процессоров 1 2 3 4 10

Время на вычисление элементов матрицы 19,060512 12,603949 8,790949 6,272236 3,200899

Время на решение СЛАУ 0,3195 0,2264 0,15966 0,11285 0,05601

Время на одну итерацию с учетом времени на передачу между процессами после итерации 0,005153226 0,003652 0,002575 0,00182 0,0009034

Общее время расчета 19,380012 12,83035 8,950609 6,385086 3,256909

С учетом времени расчета при модельных запусках при т = 4, п = 3 определено гарантированное время расчета при т = 8, п = 3 - требование при запуске на кластере комплекса СКИФ МГУ (см. выше).

Произведен модельный запуск программы при количестве процессоров, равном 100 (т = 8, п = 3), для определения времени расчета при большем значении п . Данные о времени выполнения модельного запуска программы представлены в табл. 5. Анализ результатов показывает, что резко увеличивается время выполнения программы при увеличении т:

1) за счет большего количества вычисляемых коэффициентов матрицы

2) за счет большего количества передаваемых коэффициентов матрицы и элементов векторов при решении СЛАУ (5).

С учетом времени модельного расчета при т = 8, п = 10 определено гарантированное время расчета при т = 8, п = 100.

Таблица 5

Данные о времени выполнения программы при запусках (т = 4, п = 3; т = 8, п = 3; т = 8, п = 10)

Значения констант т = 4, п = 3 т = 8, п = 3 т = 8, п = 10

Количество процессов 10 100 100

Время на вычисление элементов матрицы, с 3,200899 54,817115 26689,46108

Время на решение СЛАУ, с 0,05601 163.1978 194

Время на одну итерацию с учетом времени на передачу между процессами после итерации, с 0,000903387 0,722114 0,678322

Общее время расчета, с 3,256909 218,0149 26883,46

Общее время расчета, мин 0,054281817 3,633582 448,0577

Общее время расчета, ч 0,000904697 0,06056 7,467628

Произведен запуск программы при количестве процессоров, равном 100, т = 8, п = 100. Данные о времени выполнения программы представлены в табл. 5. Анализ результатов показывает, что увеличилось время вычисления коэффициентов матрицы за счет увеличения количества узлов интегрирования при вычислении каждого коэффициента.

По результатам работы сделаны следующие выводы:

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

2. Необходимо учитывать, что возможности любого суперкомпьютер-ного комплекса ограничены. Для поставленной задачи в основном «узким» местом является объем оперативной памяти.

3. До запуска основной программы-задачи для определения ее времени выполнения на суперкомпьютерном комплексе СКИФ МГУ необходимы запуски модельных задач.

4. Необходимо учитывать, что с увеличением количества передаваемых элементов увеличивается время на передачу элементов.

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

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

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

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

3. Описание суперкомпьютера СКИФ МГУ [Электронный ресурс]. - Режим доступа: http://parallel.ru/cluster/skif_msu.html

4. Средства программирования параллельных задач [Электронный ресурс]. - Режим доступа: http:IIwww.parallel.ru

5. MPI: A Message - Passing Interface Standart. Version i.0. - University of Tennessee. -І994. - May, 5.

6. Система управления заданиями. Запуск задач на кластере [Электронный ресурс]. -Режим доступа: http:IIwww.parallel.ru