авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

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

УДК 519.63

ПАКЕТ ПАРАЛЛЕЛЬНЫХ ПРИКЛАДНЫХ

ПРОГРАММ HELMHOLTZ3D1

Д.С. Бутюгин

В работе представлен пакет параллельных прикладных программ Helmholtz3D, кото-

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

зависимостью от времени, распространяющиеся в трехмерных областях со сложной геомет-

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

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

Ключевые слова: Уравнения Максвелла, метод конечных элементов, итерационные ал горитмы, параллельные алгоритмы.

Введение В настоящее время существует ряд коммерческих и открытых пакетов, позволяющих проводить расчеты гармонических электромагнитных полей. К ним относятся как пакеты, ориентированные на решение именно задач электромагнетизма, например ANSYS HFSS и CST Microwave Studio, так и пакеты общего назначения для решения различных крае вых задач для дифференциальных уравнений в частных производных (Partial Dierential Equations, PDE), к которым можно отнести, например, пакет FEniCS. Однако использо вание таких пакетов может быть сопряжено с определенными сложностями. Во-первых, отсутствие подходящего вычислительного инструментария в продукте может привести к невозможности адекватно описать модель и вычислить электромагнитные поля в расчетной области. Во-вторых, закрытый характер коммерческих программных продуктов может за труднить интеграцию с уже существующими у пользователей программными комплексами.

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

В работе предлагается пакет параллельных прикладных програм Helmholtz3D. Данный пакет позволяет проводить расчеты трехмерных электромагнитных полей с гармонической зависимостью от времени, распространяющиеся в трехмерных областях со сложной гео метрией. Цель разработки пакета создание средства решения задач трехмерного элек тромагнетизма, основанного на самых современных достижения в этой области. Ключевы ми особенностями данного проекта являются следующие. Решается комплексное векторное Статья рекомендована к публикации программным комитетом Международной научной конференции Па раллельные вычислительные технологии – 18 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин уравнение Гельмгольца (2), полученное непосредственно из уравнений Максвелла без ис пользования приближений. Это позволяет использовать полученные алгоритмы при рас чете широкого спектра моделей при различных физических параметрах, в частности, рас считывать электромагнитные поля в широком спектре частот. Использование различных вариационных формулировок задачи позволяет выбирать оптимальный алгоритм для про ведения вычислений в зависимости от физических параметров задачи и дает возможность получить физически корректную картину распределения электромагнитных полей. При менение символьных вычислений и техник оптимизации выражений и автогенерации кода для построения конечно-элементных (МКЭ) аппроксимаций высоких порядков (вплоть до O(h4 )) упрощает сопровождение кода и существенно снижает вероятность появления оши бок в сложных выражениях для вычисления матриц СЛАУ.

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

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



Структура данной работы следующая. В первом разделе приводится описание мате матической постановки решаемых пакетом Helmholtz3D задач. В разделе 2 приводится описание архитектуры программного комплекса. Раздел 3 содержит описание особенностей используемых технологий параллельного программирования. В разделе 4 представлены ре зультаты численных экспериментов с использованием предлагаемого пакета. В заключении обсуждаются результаты работы.

1. Математическая постановка задачи Пусть векторы напряженности электрического и магнитного полей имеют вид (Ea eit ), (Ha eit ), где t время, круговая частота, i мнимая единица, (x) действительная часть x, а Ea и Ha зависящие только от пространственных координат комплексные амплитуды. Тогда система уравнений Максвелла, описывающая электромаг нитные поля с гармонической зависимостью от времени, при отсутствии сторонних электри ческих объемных зарядов и магнитной проводимости, записывается следующим образом:

Ea = iµHa, · (r Ea ) = 0, (1) Ha = i Ea + J, · (µr Ha ) = 0.

2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d Здесь J амплитуда внешних гармонических токов, = 0 r + ie /, r = /0, µ = µ0 µr, 0 и µ 0 диэлектрическая и магнитная проницаемости вакуума, r и µr относительная диэлектрическая и магнитная проницаемости среды, а e электрическая проводимость.

Мы полагаем, что физические параметры сред не зависят от полей Ea и Ha. Тогда ли нейная система дифференциальных уравнений (1) после исключения вектора Ha сводится к уравнению Гельмгольца µ1 (2) Ea k0 r Ea = ik0 Z0 J, r где k0 = µ0 0 = /c 0, c скорость света, Z0 = µ0 /0 = µ0 c.

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

(3) n Ea = n E0, n Ha = n H0, где n внешняя нормаль к границе, а E0, H0 заданные функции координат.

Будем полагать, что расчетная область состоит из m непересекающихся подобластей m k, в каждой из которых физические параметры сред r, µr и e являются доста = k= точно гладкими. Полагаем, что внутренние границы их раздела k,k = k k являются липшицевыми и, соответственно, на них почти всюду определена нормаль nk,k, направлен ная из k в k. Выполняются условия сопряжения nk,k · (k Ea,k k Ea,k ) = 0, nk,k (Ea,k Ea,k ) = 0, (4) nk,k · (µk Ha,k µk Ha,k ) = 0, nk,k (Ha,k Ha,k ) = 0.

Предполагаем, что k0 не является максвелловским собственным числом ( не является резонансной частотой), т.е. краевая задача (2)–(4) при E0 = 0, H0 = 0 и J = 0 имеет только нулевое решение.

Вводятся стандартные соболевские пространства:

3 3 H 1 = L2 () : L2 (), H rot = L2 () L2 () :, H0 = H rot : n rot =0.

Мы полагаем, что существует E H rot такая, что n E = n E0. Введем билиней ные формы s, m : H rot H rot C:

µ1 ( u) · ( v) d, r (u · v) d, s(u, v) = m(u, v) = r и линейный функционал L : H rot C L(v) = k0 m(E, v) s(E, v) ik0 Z0 J · vd ik0 Z0 H0 · (n v) d.

Вариационная формулировка уравнения (2) с краевыми условиями (3) в форме Галер кина имеет следующий вид (см. [1, 2]): найти такое E H0, что для всех H rot rot выполнено (5) s(E, ) k0 m(E, ) = L(), 20 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин при этом исходное решение Ea восстанавливается как Ea = E + E.

Рассмотрим теперь дискретный аналог вариационной задачи. Для простоы мы рас сматриваем только тетраэдральные неструктурированные сетки. Пусть построено соответ ствующее разбиение расчетной области на непересекающиеся тетраэдральные элементы = T. В каждом из тетраэдров вводятся базисные функции, соответствующие его сте пеням свободы. Пусть Wl конечномерное пространство базисных функций порядка не выше l, конформное H rot. В работе рассматриваются иерархические базисные функции, предложенные в [3]:

Wl = W1 W2... Wl, W1 = A1, Wi = Ai Vi, где Wi инкрементальное подпространство с базисными функциями порядка i, при этом i инкрементальное подпространство со скалярными базисными функциями, конформ V ными H 1, а Ai инкрементальное подпространство с роторными базисными функциями.

Подпространства Vl,0 и Wl,0 вводятся естественным образом.

Приближенное решение E h будем искать в виде Eh = ui i.





i Конечно-элементая функция E строится таким образом, что коэффициенты разложения h по базисным функциям с ненулевым следом на принимают фиксированные значения, соответствующие первому краевому условию в (3), а остальные равны нулю.

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

0 0 0 0 (6) Mi,j = m(j, i ), Si,j = s(j, i ), fi = L(i ), где базисные функции i, j Wl,0. Тогда итоговая система принимает вид 0 (7) S k0 M u = f.

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

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

1) геометрическое и функциональное моделирование расчетной области;

2) генерация сетки в расчетной области;

3) построение МКЭ аппроксимации и генерация СЛАУ;

4) решение СЛАУ с разреженными матрицами высоких порядков;

5) постобработка решения и его визуализация.

2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d Пакет Helmholtz3D ориентирован на решение этапов 2–5. Геометрическое и функцио нальное моделирование могут быть осуществлены средствами CAD-систем. Для этого па кет предоставляет возможность импорта сеток из сторонних пакетов за счет использования подходящих конвертеров.

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

Рис. 1. Архитектура пакета Helmholtz3D В качестве языка реализации был выбран C++. Использование шаблонов (templates) в сочетании с объектно-ориентированным подходом позволяет писать максимально обобщен ные коды, не теряющие скорость работы при использовании современных компиляторов.

В качестве компилятора использовался GCC 4.4.1. Для упрощения разработки аппрокси матора и алгебраических решателей использовалась C++ библиотека линейной алгебры Eigen версии 3.0.1 [5]. В качестве основных достоинств библиотеки можно привести следу ющие: универсальность, высокая скорость работы за счет явной векторизации и оптимиза ции выражений, простота и выразительность пользовательского программного интерфейса (API), а также распространение под лицензией MPL. Помимо Eigen, использовалась высо копроизводительная библиотека Intel R Math Kernel Library (Intel R MKL [6]). Из данной библиотеки использовался прямой решатель PARDISO для разреженных СЛАУ.

2.1. Генерация сеток Задача генерации сеток является очень важной, так как сама сетка оказывает боль шое влияние на точность полученного решения. Такие сетки должны иметь адаптивные сгущения как вблизи мелких геометрических объектов расчетной области, так и вблизи особенностей решения, вызванными наличием высоко-контрастных сред. Данной проблеме посвящено большое количество работ отечественных и зарубежных авторов, однако эта про блема выходит за рамки текущей работы. Поэтому в пакете Helmholtz3D был разработан простейший генератор сеток, имеющий возможность построения квазиструктурированных сеток в областях, составленных из параллелепипедов. Однако для того, чтобы пользовате ли имели возможность использования качественных адаптивных сеток, был реализован их импорт из формата пакета NETGEN [7].

22 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин 2.2. Декомпозиция расчетной области Для решения задачи на системах с распределенной памятью требуется распределить данные между узлами системы, по возможности минимизируя объем коммуникационных данных. Основной проблемой здесь является то, что алгоритмы, генерирующие разбие ния высокого качества, сами требуют больших вычислительных ресурсов. Поэтому, как правило, приходится идти на компромисс и выбирать алгоритмы, позволяющие получить достаточно хорошее разбиение за приемлемое время.

В данной работе предлагается геометрический подход к декомпозиции области (см. по дробнее [8]), заключающийся в декомпозиции сетки расчетной области при помощи постро ения упрощенного BSP-дерева [9]. Упрощение состоит в том, что секущие плоскости пред лагается проводить ортогонально осям координат. Алгоритм рекурсивно разделяет имею щееся множество тетраэдров на две приблизительно равные части, минимизируя при этом количество тетраэдров, принадлежащих обоим подобластям. При эффективной реализации и использовании достаточно регулярных сеток можно добиться алгоритмической сложно сти метода O(NT log NT ), где NT количество тетраэдров сетки.

2.3. Построение аппроксимаций вариационных постановок задач Дискретизация вариационной постановки задачи приводит к задаче (6)–(7). Как уже было отмечено, построение СЛАУ такого вида может быть проведено на основе поэле ментной технологии с вычислением локальных матриц и векторов. Дополнительный плюс такого подхода заключается в том, что части распределенной СЛАУ могут генерироваться непосредственно на тех узлах системы, на которых они в дальнейшем будут нужны. Это позволяет избавиться от необходимости хранить глобальную СЛАУ на каком-либо из узлов, что снимает ограничения на размер решаемых систем.

Базисные функции пространств Wl могут быть выписаны с использованием кусочно линейных функций i (r) координат, равных единице в одной из вершин сетке и нулю в остальных [3]. Далее, несложно заметить, что функции i 0и 0 могут быть приведены i к следующему каноническому виду:

p p p p (8) cj 1 j,1 2 j,2 3 j,3 4 j,4 Fj (1, 2, 3, 4 ), j где cj некоторая константа, 1,..., 4 четыре не равных нулю функций i в данном тетраэдре, pj,k целочисленные степени, Fj (1, 2, 3, 4 ) константная вектор-функция, не зависящая от координат. Тогда интегралы для вычисления коэффициентов матриц S и M в каждом из тетраэдров могут быть приведены к виду p pp (9) cj ck Fj (1, 2, 3, 4 ) · Fk (1, 2, 3, 4 ) j+k i11 i22... ij+k d, j k T Такие интегралы могут быть вычислены аналитически [1]. Для автоматизации про цесса генерации локальных матриц по формуле (9) был реализован модуль, по заданному описанию базисных функций генерирующий выражения для вычисления элементов мат риц и проводящий символьные оптимизации получаемых выражений. Схема программы представлена на рис. 2.

2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d EBNF Canonic Code parser transformation generator Optimizer Рис. 2. Схема работы генератора выражений Правила описания базисных функций задают некоторую формальную грамматику и могут быть записаны в расширенной форме Бэкуса–Наура (РБНФ, или EBNF). Набор ба зисных функций читается парсером РБНФ (EBNF parser) и преобразуется в промежуточ ный вид, соответствующий дереву разбора выражений. Поскольку не любое выражение, записанное в предложенной грамматике является корректным, то в парсер дополнительно встроен валидатор, который проверяет выражение на совместимость.

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

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

2.4. Распределенные итерационные решатели После декомпозиции задачи на подобласти и построения блоков матрицы на соответ ствующих узлах кластера, мы имеем распределенную расширенную СЛАУ u A = f.

Для решения таких СЛАУ широко используются методы Шварца. Одна итерация аддитив ного метода Шварца в матричном виде может быть записана как (см. [10]) Rp A1 Rp (f Aui ), T ui+1 = ui + p p где Rp оператор сужения на подобласть p, p = 1,..., P ;

здесь и далее подчеркивание для ui мы опускаем. Матрицу Ap для подобласти p можно представить в виде Ap = Rp ARp. T Проблема состоит в том, что для уравнения Гельмгольца метод Шварца может не схо диться. Однако можно заметить, что решение системы методом Шварца есть ни что иное, как решение методом простой итерации предобусловленной системы u B 1 A = B 1 f, B 1 = Rp A1 Rp.

T (10) p p Такую предобусловленную систему можно решать итерационными методами в подпро странствах Крылова. Однако метод Шварца имеет недостаток, заключающийся в том, что 24 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин количество итераций, необходимое для достижения заданной точности, растет с увеличе нием числа подобластей [10]. Для преодоления этого недостатка можно воспользоваться методом коррекции на грубой сетке [11] (Coarse Grid Correction, CGC).

Для этого введем операторы грубосеточного сужения Rc, определенный как оператор сужения на базис подпространства W1, а также оператор продолжения Pc = Rc. Тогда T итерация метода аддитивного Шварца с грубосеточной коррекцией записывается как ui+1/2 = ui + Pc A1 Rc (f Aui ), ui+1 = ui+1/2 + p Rp p p i+1/2 ), T A1 R (f Au где A1 = Rc APc матрица, полученная при использовании базисных функций порядка 1.

Мы будем решать предобусловленную систему вида (10), где предобуславливатель B при наличии грубосеточной коррекции имеет вид Rp A1 Rp I APc A1 Rc + Pc A1 Rc.

B 1 = T (11) p 1 p Поскольку метод коррекции на грубой сетке не очень чувствителен к точности аппрок симации A1 [12], для приближенного решения систем вида A1 v = q можно использовать итерационные методы в подпространствах следов, см. например [13]. Однако приближенное обращение матрицы A1 приводит к тому, что предобуславливатель B 1 становится пере менным. В этом случае хорошим выбором для итерационного решения задач с таким пре добуславливателем оказывается метод Flexible GMRES [10], допускающий использование переменных предобуславливателей.

Для решения задач в подобластях Ap x = b можно использовать какой-нибудь пря мой решатель для разреженных систем, например решатель PARDISO из библиотеки Intel MKL [6]. Проблема состоит в том, что время, требуемое PARDISO для разложения матриц, растет практически как O(N 3 /P 3 ), поэтому такой метод подходит только для решения не очень больших систем на большом числе процессоров. Однако действие обратных матриц A1 можно также вычислять приближенно, например методом FGMRES с предобуславли p вателем AMG, основанным на использовании иерархических базисов, предложенным в [13].

2.5. Методы визуализации электромагнитных полей Визуализация электромагнитных полей является неотъемлемой частью численного экс перимента, поскольку позволяет проанализировать характерное распределение полей в рас четной области. В связи с этим, в рамках пакета Helmholtz3D был реализован модуль ви зуализации электромагнитных полей, основанный на технологии OpenGL.

Визуализация электромагнитных полей осуществляется следующим образом. Визуали затор предоставляет возможность визуализировать одно из полей – электрическое либо магнитное, и позволяет переключаться между этими двумя режимами. Поля отрисовыва ются в каждом из тетраэдров в виде стрелки, сонаправленной с полем в данной точке, и имеющей длину, пропорциональную амплитуде поля в текущий момент времени. Для боль шей наглядности амплитуда поля кодируется цветом – с линейным переходом от желтого (максимально возможная амплитуда для текущего поля), до синего (близкая к нулю ам плитуда). Отметим, что поскольку амплитуда поля представляется в виде двух векторов – действительной и мнимой компоненты, то в заданный момент времени t реальная ампли 2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d туда поля задается формулой F (t) = Fre cos(t) + Fim sin(t), где Fre и Fim векторы амплитуд действительной и мнимой части для соответствующего поля (электрического либо магнитного). Для примера на рис. 3 приведена визуализация электрического поля в волноводе с однородной средой. На рисунке хорошо видна структура электрического поля.

Рис. 3. Пример визуализации электрического поля в волноводе 3. Технологии параллельного программирования и оптимизации для систем с распределенной и общей памятью Параллельные распределенные алгоритмы в пакете Helmholtz3D построены на основе широко распространенного стандарта интерфейса обмена данными MPI (Message Passing Interface). При этом используемые в пакете итерационные решатели с предобуславливателя ми типа аддитивного Шварца естественным образом ложатся на архитектуру современных вычислительных систем. Действительно, итерационное решение СЛАУ методом FGMRES требует 3 следующих типа операций.

К первому типу относятся векторно-векторные операции. Для имеющейся декомпози ции задачи на подобласти распараллеливание таких операций проводится естественным образом, за исключением операций вычисления скалярного произведения и нормы, для ко торых необходима точка синхронизации и обмен данными. Однако в интерфейсе MPI имеет ся требуемая функция MPI_Allreduce, реализация которой оптимизируется поставщиком библиотеки MPI.

Для умножения матрицы на вектор решатель строит в каждой подобласти списки уз лов, являющихся граничными для соседних подобластей. В дальнейшем решатель исполь зует собранную информацию для пересылки частей вектора между процессами. Ключевым моментом в таком подходе является уменьшение объема пересылаемых данных в сосед ние подобласти требуется пересылать только граничные коэффициенты векторов вместо того, чтобы пересылать вектор целиком. В итоге, общий объем пересылаемых данных в трехмерных задачах электромагнетизма для одного умножения матрицы на вектор будет пропорционален O(N 2/3 log P ) вместо O(N ).

Применение предобуславливателей, основанных на методе Шварца, требует на каждой итерации внешнего алгоритма решения подзадач в подобластях, и, согласно (10), не требу ют коммуникаций. Распараллеливание решателя в подобластях проводилось с использова 26 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин нием средств OpenMP. Были реализованы параллельные процедуры умножения матрицы на вектор, а также параллельное решение разреженных треугольных систем для операто ра сглаживания SSOR (Symmetric Successive Over-Relaxation, [10]) в предобуславливателе AMG. Распараллеливание треугольного решателя проводилось при помощи метода плани рования вычислений по уровням: при решении системы вида U x = y для каждой вершины вычисляется ее уровень li :

(12) li = max {lj + 1 : Ui,j = 0, j i}.

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

Можно отметить ряд полезных оптимизаций решателя в подобластях. Для начала необ ходимо отметить его особенность: такое итерационное решение СЛАУ в подобластях оказы вается bandwidth-limited, то есть существенно ограничено пропускной способностью шины данных системы, поскольку за одну итерацию алгоритма, фактически, требуется несколько раз прочитать всю матрицу системы (в том числе для решения двух треугольных систем), и несколько раз прочитать и записать различные векторы. Таким образом, реализован ные оптимизации можно разделить на две категории: оптимизации доступа к памяти за счет переупорядочивания матрицы и оптимизации, направленные на повышение произво дительности алгоритма на NUMA-системах, заключающиеся в правильной инициализации областей памяти. Подробнее об этих оптимизациях см. [14].

Отметим, что помимо матрично-векторных операций, решатель в подобластях исполь зует прямой метод PARDISO из Intel R MKL для самой грубой системы в предобуслав ливателе AMG. Дополнительное распараллеливание и оптимизация данной части решате ля не требуется: PARDISO поддерживает OpenMP параллелизм как на этапе построения LU -разложения матрицы, так и на этапе поиска решения, а также использует высокоэф фективные функции BLAS и LAPACK.

4. Численные эксперименты Предлагаемые в работе алгоритмы были апробированы на серии методических и прак тических задач электромагнетизма.

Первой из задач является вычисление электромагнитного поля в параллелепипедаль ном волноводе [0, a] [0, b] [0, c], a = 72 мм, b = 34 мм и c = 200 мм, см. рис. 4.

Рис. 4. Расчетная область для модельной задачи Физические параметры среды: µr = 1, r = 1 0.1i, = 6 · 109 Гц. При z = 200 мм поставлено краевое условие (3) c E0 = ey sin(x/a), на остальной части границы краевое 2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d условие идеального проводника E0 = 0. Известен аналитический вид решения:

x sin z (13) E = ey sin.

a sin c Сетка в данном случае была сгенерирована средствами пакета Helmholtz3D.

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

Параметры µr = 1, r = 1 во всей области. Проводимость среды = 0.1 См/м, в среде имеется слой с = 0.05 См/м. Координаты слоя по оси Oz: от 0.425 м до 0.275 м. В среде пробурена вертикальная цилиндрическая скважина радиуса 0.108 м с центром в на чале координат в плоскости Oxy. В скважине имеется цилиндрическая каверна внешнего радиуса 0.118 м и положением по оси Oz от 0.0725 до 0.0725 м. Проводимость скважины и каверны = 5 См/м. В скважину вставлен полый цилиндрический прибор радиуса 0. м, смещенный по оси x относительно центра скважины на 0.064, заполненный воздухом ( = 0). В приборе имеется одна генераторная и две приемные петли радиуса 0.0365 м при z = 0.5 м и при z = 0.0 и z = 0.1 м соответственно. По генераторной петле течет ток 1 А с частотой 14 МГц. Искомой является разность фаз ЭДС в приемниках.

Генерация неравномерной сетки проводилась при помощи утилиты NETGEN [7]. Коли чество тетраэдров сетки NT = 490282, порядок СЛАУ без учета пересечений N0 = 8946864.

Тестирование проводились на двух системах: OpenMP решатель для систем с общей памятью (соответствует P = 1) на системе Intel R Xeon X7560, 2.27 ГГц, 4 8 ядер, 64 ГБ памяти, гибридный решатель на узлах HP BL2х220 G6 кластера НКС-30Т. Каждый узел кластера Intel Xeon Е5540, 2.53 ГГц, 2 4 ядер, 16 ГБ памяти.

Для аппроксимации использовались базисные функции W3 третьего порядка.

В качестве критерия останова итерационного процесса было выбрано уменьшение нор мы невязки в раз: |ri | |f |. Для внутренних итераций FGMRES с предобуславливателем AMG A = 0.2, для грубосеточной коррекции в подпространстве следов c = 0.1. Решение задач в подобластях при обращении Ap проводилось с p = 0.01. Внешние итерации прово дились с точностью = 107 за исключением сетки 116 55 320 для первой задачи, для которой использовалось = 109 для получения E h нужной точности.

В табл. 1–3 приведены результаты тестирования решателей на модельной задаче и за даче электромагнитного каротажа: N0 порядок системы без учета пересечений, E мак Используется более одного MPI процесса на узел (всего 64 узла).

28 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин Таблица Масштабируемость гибридного кластерного решателя на модельной задаче Сетка, N0 Количество подобластей (MPI процессов) 1281 2561 и E 16 32 705312 1168248 1777626 2599050 3928914 58 28 160 Nb 8 9 10 10 11 n 28519914 266 385 382 468 577 nc 24 27 29 29 32 nA 1.2e3 5.2e2 2.3e2 1.1e2 8.0e1 5.7e 7,1 · 107 t 13035042 116 55 320 Nb 13 n 225335973 N/A N/A N/A N/A 1102 nc 40 nA 1.2e3 9.2e 9,0 · 108 t симальная относительная ошибка электрического поля, Nb размерность пространства следов для A1, n количество внешних итераций распределенного решателя, n3 и n общее количество внешних и внутренних (в предобуславливателе AMG) итераций OpenMP решателя, nc общее количество итераций грубосеточной коррекции, nA максимальное количество внешних итераций метода FGMRES в подобластях, t время решения СЛАУ в секундах в экспоненциальном формате a.ben = a.b · 10n. Результат решения задачи электро магнитного каротажа показал хорошее соответствие с результатом группы НГТУ, которая решала эту задачу другим методом методом выделения поля источника, при этом иско мая разность фаз ЭДС в приемниках совпала с ошибкой 1.8 %.

Таблица Масштабируемость OpenMP решателя на задаче электромагнитного каротажа Количество OpenMP потоков 1 2 4 8 16 14 14 14 14 14 n 33 33 33 33 33 n 7.5e3 4.1e3 2.2e3 1.3e3 8.8e2 6.1e t Из таблиц видно, что предобуславливатели AMG и CGC успешно уменьшают коли чество внешних итераций до приемлемого значения. Количество итераций грубосеточной коррекции растет, но гораздо медленнее, чем линейно, за счет трехмерной декомпозиции расчетной области. Тем не менее, при дальнейшем росте размера СЛАУ может потребовать ся геометрическое огрубление задачи для A1, например [15], в результате чего мы получим многоуровневую схему грубосеточной коррекции.

Используется более одного MPI процесса на узел (всего 64 узла).

2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d Таблица Масштабируемость MPI решателя на задаче электромагнитного каротажа Количество подобластей (MPI процессов) 1281 2561 4 8 16 32 102918 236787 477672 803112 1346418 2038260 2873652 Nb 14 19 19 36 35 36 38 n 86 131 148 310 310 341 383 nc 49 61 61 111 109 113 123 nA 8.8e2 4.9e2 2.3e2 1.9e2 9.2e1 6.2e1 5.5e1 4.9e t Заключение В работе представлен пакет параллельных прикладных програм Helmholtz3D. Данный пакет позволяет проводить расчеты трехмерных электромагнитных полей с гармонической зависимостью от времени. Предложенные алгоритмы позволяют рассчитывать электромаг нитные поля в с областях со сложной разномасштабной геометрией и высококонтрастными средами с высокой точностью, продемонстрированной на методических и практических за дачах. Параллельные алгоритмы решения СЛАУ для систем с общей и распределенной памятью позволяют решать сверхбольшие задачи до нескольких сотен миллионов неиз вестных за время порядка 10–15 минут.

Работа поддержана грантами РФФИ №11-01-00205 и ОМН РАН №1.3.3-4.

Литература 1. Соловейчик, Ю.Г. Метод конечных элементов для решения скалярных и векторных задач / Ю.Г. Соловейчик, М.Э. Рояк, М.Г. Персова. Новосибирск: Изд-во НГТУ, 2007.

2. Monk, P. Finite Element Methods for Maxwell’s Equations / P. Monk. Oxford University Press, 2003.

3. Ingelstrm, P. A new set of H(curl)-conforming hierarchical basis functions for tetrahedral o meshes / P. Ingelstrm // IEEE Transactions on Microwave Theory and Techniques. 2006.

o Vol. 54, № 1. P. 160–114.

4. Ильин, В.П. Методы и технологии конечных элементов / В.П. Ильин. Новосибирск:

Изд-во ИВМиМГ СО РАН, 2007.

5. Eigen. URL: http://eigen.tuxfamily.org (дата обращения 12.12.2012).

6. Intel (R) Math Kernel Library from Intel. URL: http://software.intel.com/en-us/ intel-mkl (дата обращения 12.02.2013).

7. Schberl, J. NETGEN An advancing front 2D/3D-mesh generator based on

Abstract

rules o / J. Schberl // Computing and Visualization in Science. 1997. Vol. 1, № 1. P. 41–52.

o 8. Бутюгин, Д.С. Алгоритмы решения СЛАУ на системах с распределенной памятью в применении к задачам электромагнетизма / Д.С. Бутюгин // Вестник ЮУрГУ. Серия Вычислительная математика и информатика. 2012. № 46(305). С. 5–18.

30 Вестник ЮУрГУ. Серия Вычислительная математика и информатика Д.С. Бутюгин 9. Fuchs, H. On visible surface generation by a priori tree structures / H. Fuchs, Z.M. Kedem, B.F. Naylor // ACM Computer Graphics. 1980. Vol. 14, № 3. P. 124–133.

10. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition / Y. Saad. SIAM, 2003.

11. Bramble, J. The construction of preconditioners for elliptic problems by substructuring. I.

/ J. Bramble, J. Pasciak, A. Schatz // Mathematics of Computation. 1986. Vol. 47, № 175. P. 103–134.

12. Nabben, R. A comparison of deation and coarse grid correction applied to porous media ow / R. Nabben, C. Vuik // SIAM Journal on Numerical Analysis. 2004. Vol. 42, № 4.

P. 1631–1647.

13. Butyugin, D.S. Ecient iterative solvers for time-harmonic Maxwell equations using domain decomposition and algebraic multigrid / D.S. Butyugin // Journal of Computational Science.

2012. Vol. 3, № 6. P. 480–485.

14. Бутюгин, Д.С. Параллельный предобуславливатель SSOR для решения задач электро магнетизма в частотной области / Д.С. Бутюгин // Вычислительные методы и про граммирование. 2011. Т. 12, № 1. С. 110–117.

15. Toward an h-independent algebraic multigrid method for Maxwell’s equations / J. Hu, R.

Tuminaro, P. Bochev et al. // SIAM Journal on Scientic Computing. 2005. Vol. 27, № 5. P. 1669–1688.

Дмитрий Сергеевич Бутюгин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, аспирант, Новосибирский государствен ный университет, dm.butyugin@gmail.com.

PARALLEL APPLICATION PACKAGE HELMHOLTZ3D D.S. Butyugin, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation) The paper present parallel application package Helmholtz3D, which allows to compute three-dimensional time-harmonic electromangetic elds propagating in domains with complicated geometry. In order to solve systems of linear algebraic equations (SLAEs) with complex ill conditioned non-hermitian matrices arising in the nite element discretizations, the package uses state-of-the-art iterative methods in Krylov subspace combined with original parallel preconditioners. Package approbation was performed on the series of model and real-life problems of electromangetic eld computation in the microwave devices and well logging problems.

Keywords: Maxwell’s equations, nite element method, iterative algorithms, parallel algorithms.

References 1. Soloveychick Y.G., Royak M.E., Persova M.G. Metod konechnykh elementov dlya resheniya skalarnykh i vectornykh zadach [Finite Element Method for the Solution of Scalar and Vector Problems]. Novosibirsk, NSTU Publ., 2007.

2. Monk P. Finite Element Methods for Maxwell’s Equations. Oxford University Press, 2003.

2013, т. 2, № Пакет параллельных прикладных программ Helmholtz3d 3. Ingelstrm P. A new set of H(curl)-conforming hierarchical basis functions for tetrahedral o meshes. IEEE Transactions on Microwave Theory and Techniques. 2006. Vol. 54, No. 1.

P. 160–114.

4. Il’in V.P. Metody i tekhnologii konechnykh elementov [Methods and technologies of nite elements]. Novosibirsk, ICM&MG SBRAS Publ., 2007.

5. Eigen. URL: http://eigen.tuxfamily.org (accessed 12 December 2012).

6. Intel (R) Math Kernel Library from Intel. URL: http://software.intel.com/en-us/ intel-mkl (accessed 12 Februrary 2013).

7. Schberl J. NETGEN An advancing front 2D/3D-mesh generator based on abstract rules.

o Computing and Visualization in Science. 1997. Vol. 1, No. 1. P. 41–52.

8. Butyugin D.S. Algoritmy resheniya SLAU na sistemakh c raspredelennoy pamatyu v primemenii k zadacham electromagnetisma [Algorithms of SLAEs solution on the systems with distributed memory applied to the problems of electromagnetism]. Vestnik YUURGU.

Seriya “Vychislitelnaya matematika i informatika” [Bulletin of South Ural State University.

Seriers: Computational Mathematics and Software Engineering], 2012. No. 46(305). P. 5–18.

9. Fuchs H., Kedem Z.M., Naylor B.F. On visible surface generation by a priori tree structures ACM Computer Graphics. 1980. Vol. 14, No. 3. P. 124–133.

10. Saad Y. Iterative Methods for Sparse Linear Systems, Second Edition. SIAM, 2003.

11. Bramble J., Pasciak J., Schatz A. The construction of preconditioners for elliptic problems by substructuring. I. Mathematics of Computation. 1986. Vol. 47, No. 175. P. 103–134.

12. Nabben R., Vuik C. A comparison of deation and coarse grid correction applied to porous media ow. SIAM Journal on Numerical Analysis. 2004. Vol. 42, No. 4. P. 1631–1647.

13. Butyugin D.S. Ecient iterative solvers for time-harmonic Maxwell equations using domain decomposition and algebraic multigrid. Journal of Computational Science. 2012. Vol. 3, No. 6.

P. 480–485.

14. Butyugin D.S. Parallelniy predobuslavlivatel’ SSOR dlya resheniya zadach electromagnetisma v chastotnoy oblasty [Parallel SSOR preconditioner for the solution of electromagnetism problems in the frequency domain]. Vychislitelnye metody i programmirovanie [Computational methods and programming]. 2011. Vol. 12, No. 1.

P. 110–117.

15. Hu J., Tuminaro R., Bochev P., Garasi C., Robinson C. Toward an h-independent algebraic multigrid method for Maxwell’s equations SIAM Journal on Scientic Computing. 2005.

Vol. 27, No. 5. P. 1669–1688.

Поступила в редакцию 15 февраля 2013 г.

32 Вестник ЮУрГУ. Серия Вычислительная математика и информатика

 

Похожие работы:





 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.