CC BY
115
ПАКЕТ ПАРАЛЛЕЛЬНЫХ ПРИКЛАДНЫХ ПРОГРАММ HELMHOLTZ3D1 Д.С. Бутюгин
В работе представлен пакет параллельных прикладных программ Helmholtz3D, который позволяет проводить расчеты трехмерных электромагнитных полей с гармонической зависимостью от времени, распространяющиеся в трехмерных областях со сложной геометрией. Для решения возникающих в результате аппроксимаций систем линейных алгебраических уравнений (СЛАУ) с комплексными плохообусловленными неэрмитовыми матрицами используются современные итерационные методы решения СЛАУ в подпространствах Крылова совместно с оригинальными параллельными предобуславливателями. Апробация пакета проведена на серии методических и практических задач расчета электромагнитных полей для волновых устройств и задач электромагнитного каротажа.
Ключевые слова: Уравнения Максвелла, метод конечных элементов, итерационные алгоритмы, параллельные алгоритмы.
Введение
В настоящее время существует ряд коммерческих и открытых пакетов, позволяющих проводить расчеты гармонических электромагнитных полей. К ним относятся как пакеты, ориентированные на решение именно задач электромагнетизма, например ANSYS HFSS и CST Microwave Studio, так и пакеты общего назначения для решения различных краевых задач для дифференциальных уравнений в частных производных (Partial Differential Equations, PDE), к которым можно отнести, например, пакет FEniCS. Однако использование таких пакетов может быть сопряжено с определенными сложностями. Во-первых, отсутствие подходящего вычислительного инструментария в продукте может привести к невозможности адекватно описать модель и вычислить электромагнитные поля в расчетной области. Во-вторых, закрытый характер коммерческих программных продуктов может затруднить интеграцию с уже существующими у пользователей программными комплексами. В третьих, далеко не всегда имеется поддержка вычислений на удаленных вычислительных системах, имеющихся в распоряжении институтов и коммерческих компаний. И наконец, не стоит сбрасывать со счетов высокую стоимость лицензии, которая может достигать десятков тысяч долларов в год. В открытых же пакетах, зачастую, реализованы не слишком эффективные алгоритмы, а имеющаяся по ним документация ограничена. В связи с этим, актуальной является разработка пакетов, основанных на самых современных достижениях в области конечно-элементного моделирования электромагнитных полей и поддерживающих расчеты на параллельных системах с общей и распределенной памятью с тысячами и десятками тысяч вычислительных узлов.
В работе предлагается пакет параллельных прикладных програм Helmholtz3D. Данный пакет позволяет проводить расчеты трехмерных электромагнитных полей с гармонической зависимостью от времени, распространяющиеся в трехмерных областях со сложной геометрией. Цель разработки пакета — создание средства решения задач трехмерного электромагнетизма, основанного на самых современных достижения в этой области. Ключевыми особенностями данного проекта являются следующие. Решается комплексное векторное
1 Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии - 2013»
уравнение Гельмгольца (2), полученное непосредственно из уравнений Максвелла без использования приближений. Это позволяет использовать полученные алгоритмы при расчете широкого спектра моделей при различных физических параметрах, в частности, рассчитывать электромагнитные поля в широком спектре частот. Использование различных вариационных формулировок задачи позволяет выбирать оптимальный алгоритм для проведения вычислений в зависимости от физических параметров задачи и дает возможность получить физически корректную картину распределения электромагнитных полей. Применение символьных вычислений и техник оптимизации выражений и автогенерации кода для построения конечно-элементных (МКЭ) аппроксимаций высоких порядков (вплоть до 0(Л,4)) упрощает сопровождение кода и существенно снижает вероятность появления ошибок в сложных выражениях для вычисления матриц СЛАУ.
Для решения возникающих в результате аппроксимаций систем линейных алгебраических уравнений (СЛАУ) с комплексными плохообусловленными неэрмитовыми матрицами используются современные итерационные методы решения СЛАУ в подпространствах Крылова совместно с оригинальными параллельными предобуславливателями. Предобуславли-ватели используют методы аддитивного Шварца и экономичной геометрической декомпозиции области, а также алгебраические мультисеточные методы на основе иерархических базисов, в том числе методы алгебраической грубосеточной коррекции. Применение таких алгоритмов позволяет добиться высокого уровня производительности решателей на системах с общей и распределенной памятью и хорошей масштабируемости на сотнях вычислительных узлов. Алгоритмы пакета были апробированы на серии методических и практических задач расчета электромагнитных полей для волновых устройств и задач электромагнитного каротажа, продемонстрировав высокую эффективность, в частности, решение получаемых СЛАУ с сотнями миллионов неизвестных может занимать около 10-15 минут. Высокая точность полученных конечно-элементных решений была показана путем сравнения с имеющимися аналитическими решениями либо результатами расчетов других групп (в том числе и полученных другими методами).
Структура данной работы следующая. В первом разделе приводится описание математической постановки решаемых пакетом Не1тИо^3Б задач. В разделе 2 приводится описание архитектуры программного комплекса. Раздел 3 содержит описание особенностей используемых технологий параллельного программирования. В разделе 4 представлены результаты численных экспериментов с использованием предлагаемого пакета. В заключении обсуждаются результаты работы.
1. Математическая постановка задачи
Пусть векторы напряженности электрического и магнитного полей имеют вид
^Я(Еае-гшЬ), ^&(Нае-гшЬ), где £ — время, ш — круговая частота, г — мнимая единица, К(ж) —
действительная часть ж, а Еа и На — зависящие только от пространственных координат комплексные амплитуды. Тогда система уравнений Максвелла, описывающая электромагнитные поля с гармонической зависимостью от времени, при отсутствии сторонних электрических объемных зарядов и магнитной проводимости, записывается следующим образом:
^ X Еа — гш^Наз V ■ (8ГЕа) — 0; (1)
V X На — —1шеЕа + , V ■ (^г На) — 0.
Здесь 3 — амплитуда внешних гармонических токов, е — е0ег + гае/ш, ег — е/е0, ц — ц0цг, ео и цо — диэлектрическая и магнитная проницаемости вакуума, ег и цг — относительная диэлектрическая и магнитная проницаемости среды, а ае — электрическая проводимость.
Мы полагаем, что физические параметры сред не зависят от полей Еа и На. Тогда линейная система дифференциальных уравнений (1) после исключения вектора На сводится к уравнению Гельмгольца
V х (д V х E^j — hferEa = ikoZoJ, (2)
где к0 — — ш/с > 0, с — скорость света, Z0 — л/ц0/е0 — ц0с.
Решение рассматриваемых уравнений будем искать в ограниченной односвязной области О с липшицевой границей дО — ГуХ, на каждой из частей которой поставлено одно из граничных условий:
ПХ Ea
= П х E0, П х Ha
г
= П х H0, (3)
s
где П — внешняя нормаль к границе, а Eo, Hо — заданные функции координат.
Будем полагать, что расчетная область Q состоит из m непересекающихся подобластей
m
Q = |^J Qk, в каждой из которых физические параметры сред ¿r, ßr и ae являются доста-k=i
точно гладкими. Полагаем, что внутренние границы их раздела T^k' = QkP|Qк являются липшицевыми и, соответственно, на них почти всюду определена нормаль nk,k', направленная из Qk в Qk'. Выполняются условия сопряжения
Пк,к' " (&kEa,k &k' Ea,k') 0, Пк,к' х (Ea,k Ea,k') 0, (4)
Пк,к' ' (ßkHa,k ßk' Ha,k') = 0, Tlk,k' х (Ha,k Ha,k') = 0.
Предполагаем, что k2 не является максвелловским собственным числом (ш не является
резонансной частотой), т.е. краевая задача (2)-(4) при Е0 = 0, Н0 = 0 и J = 0 имеет только
нулевое решение.
Вводятся стандартные соболевские пространства:
H1 = {^ е L2(Q) : Vp е [L2(Q)]3} , Hrot = {ф е [L2(Q)]3 : V х ф е [L2(Q)]3} ,
H0ot = {ф e Hrot : П х ф|г = 0} .
Мы полагаем, что существует Eг е Hrot такая, что П х Ег|г = П х Eo. Введем билинейные формы s,m : Hrot х Hrot ^ C:
s(u,v)= ß-1 (V х и) ■ ^хь)(1й, m(u,v) = ¿r (и ■v)dQ,
Jo Jo
и линейный функционал L : Hrot ^ C
L(v) = k02m(Eг, v) — s(Eг, v) — ik0Z0 J ■ vdQ — ik0Z0 H0 ■ (П х v) dQ.
Jo J s
Вариационная формулировка уравнения (2) с краевыми условиями (3) в форме Галер-кина имеет следующий вид (см. [1, 2]): найти такое E е H0ot, что для всех ф е H0ot выполнено
s(E,'ip) — k0°m(E,'H) = L(rf), (5)
при этом исходное решение Еа восстанавливается как Еа — Е + Ег.
Рассмотрим теперь дискретный аналог вариационной задачи. Для простоы мы рассматриваем только тетраэдральные неструктурированные сетки. Пусть построено соответствующее разбиение расчетной области О на непересекающиеся тетраэдральные элементы
О — и От. В каждом из тетраэдров вводятся базисные функции, соответствующие его степеням свободы. Пусть Щг — конечномерное пространство базисных функций порядка не выше I, конформное НГ°Ч В работе рассматриваются иерархические базисные функции, предложенные в [3]:
Щ — У1 ® УУ2 ® ... ®УЩ, Щ — М, Щ — Лг ФУУг,
где УУг — инкрементальное подпространство с базисными функциями порядка г, при этом
Vг — инкрементальное подпространство со скалярными базисными функциями, конформными Н1, а Лг — инкрементальное подпространство с роторными базисными функциями. Подпространства VI,0 и Щг,0 вводятся естественным образом.
Приближенное решение Ен будем искать в виде
Ен — ^ иг$.
г
Конечно-элементая функция Ер строится таким образом, что коэффициенты разложения по базисным функциям с ненулевым следом на Г принимают фиксированные значения, соответствующие первому краевому условию в (3), а остальные равны нулю.
Вводятся матрицы соответствующих билинейных форм и вектор правой части:
— т{'$,'$), — ¡г — Ь($), (6)
где базисные функции Н0,Н° € Уг,0. Тогда итоговая система принимает вид
[5 - к£м] и — (7)
Для вычисления элементов матрицы и вектора правой части можно воспользовать-
ся поэлементной технологией сборки [4], заменив интегрирование по расчетной области О суммой интегралов по каждому из тетраэдров и вычислением в каждом из тетраэдров локальных матриц и векторов.
2. Методология построения пакета параллельных прикладных программ Helmholtz3D
Задача моделирования гармонических электромагнитных полей может быть разбита на следующие функциональные этапы:
1) геометрическое и функциональное моделирование расчетной области;
2) генерация сетки в расчетной области;
3) построение МКЭ аппроксимации и генерация СЛАУ;
4) решение СЛАУ с разреженными матрицами высоких порядков;
5) постобработка решения и его визуализация.
Пакет Не1тЬоН^3Б ориентирован на решение этапов 2-5. Геометрическое и функциональное моделирование могут быть осуществлены средствами СЛБ-систем. Для этого пакет предоставляет возможность импорта сеток из сторонних пакетов за счет использования подходящих конвертеров.
В соответствии с функциональным назначением пакета была разработана его архитектура. Схематично она представлена на рис. 1.
Рис. 1. Архитектура пакета Helmholtz3D
В качестве языка реализации был выбран C++. Использование шаблонов (templates) в сочетании с объектно-ориентированным подходом позволяет писать максимально обобщенные коды, не теряющие скорость работы при использовании современных компиляторов. В качестве компилятора использовался GCC 4.4.1. Для упрощения разработки аппрокси-матора и алгебраических решателей использовалась C++ библиотека линейной алгебры Eigen версии 3.0.1 [5]. В качестве основных достоинств библиотеки можно привести следующие: универсальность, высокая скорость работы за счет явной векторизации и оптимизации выражений, простота и выразительность пользовательского программного интерфейса (API), а также распространение под лицензией MPL. Помимо Eigen, использовалась высокопроизводительная библиотека Intel® Math Kernel Library (Intel® MKL [6]). Из данной библиотеки использовался прямой решатель PARDISO для разреженных СЛАУ.
2.1. Генерация сеток
Задача генерации сеток является очень важной, так как сама сетка оказывает большое влияние на точность полученного решения. Такие сетки должны иметь адаптивные сгущения как вблизи мелких геометрических объектов расчетной области, так и вблизи особенностей решения, вызванными наличием высоко-контрастных сред. Данной проблеме посвящено большое количество работ отечественных и зарубежных авторов, однако эта проблема выходит за рамки текущей работы. Поэтому в пакете Не1тЬоН^3Б был разработан простейший генератор сеток, имеющий возможность построения квазиструктурированных сеток в областях, составленных из параллелепипедов. Однако для того, чтобы пользователи имели возможность использования качественных адаптивных сеток, был реализован их импорт из формата пакета NETGEN [7].
2.2. Декомпозиция расчетной области
Для решения задачи на системах с распределенной памятью требуется распределить данные между узлами системы, по возможности минимизируя объем коммуникационных данных. Основной проблемой здесь является то, что алгоритмы, генерирующие разбиения высокого качества, сами требуют больших вычислительных ресурсов. Поэтому, как правило, приходится идти на компромисс и выбирать алгоритмы, позволяющие получить достаточно хорошее разбиение за приемлемое время.
В данной работе предлагается геометрический подход к декомпозиции области (см. подробнее [8]), заключающийся в декомпозиции сетки расчетной области при помощи построения упрощенного BSP-дерева [9]. Упрощение состоит в том, что секущие плоскости предлагается проводить ортогонально осям координат. Алгоритм рекурсивно разделяет имеющееся множество тетраэдров на две приблизительно равные части, минимизируя при этом количество тетраэдров, принадлежащих обоим подобластям. При эффективной реализации и использовании достаточно регулярных сеток можно добиться алгоритмической сложности метода 0(Nt log Nt), где Nt — количество тетраэдров сетки.
2.3. Построение аппроксимаций вариационных постановок задач
Дискретизация вариационной постановки задачи приводит к задаче (6)—(7). Как уже было отмечено, построение СЛАУ такого вида может быть проведено на основе поэлементной технологии с вычислением локальных матриц и векторов. Дополнительный плюс такого подхода заключается в том, что части распределенной СЛАУ могут генерироваться непосредственно на тех узлах системы, на которых они в дальнейшем будут нужны. Это позволяет избавиться от необходимости хранить глобальную СЛАУ на каком-либо из узлов, что снимает ограничения на размер решаемых систем.
Базисные функции пространств Wi могут быть выписаны с использованием кусочнолинейных функций ^i(r) координат, равных единице в одной из вершин сетке и нулю — в остальных [3]. Далее, несложно заметить, что функции ф0 и V х ф0 могут быть приведены к следующему каноническому виду:
где сз — некоторая константа, £і,... ,£4 — четыре не равных нулю функций £ в данном тетраэдре, — целочисленные степени, (£1, £2, £з, £4) — константная вектор-функция,
не зависящая от координат. Тогда интегралы для вычисления коэффициентов матриц 5 и М в каждом из тетраэдров могут быть приведены к виду
Такие интегралы могут быть вычислены аналитически [1]. Для автоматизации процесса генерации локальных матриц по формуле (9) был реализован модуль, по заданному описанию базисных функций генерирующий выражения для вычисления элементов матриц и проводящий символьные оптимизации получаемых выражений. Схема программы представлена на рис. 2.
Е Е cj(F (6,6,6,Є4) ■ Fk (6,6,6,Є4)) / с с... j dn, (9)
Рис. 2. Схема работы генератора выражений
Правила описания базисных функций задают некоторую формальную грамматику и могут быть записаны в расширенной форме Бэкуса-Наура (РБНФ, или EBNF). Набор базисных функций читается парсером РБНФ (EBNF parser) и преобразуется в промежуточный вид, соответствующий дереву разбора выражений. Поскольку не любое выражение, записанное в предложенной грамматике является корректным, то в парсер дополнительно встроен валидатор, который проверяет выражение на совместимость.
Проводившиеся символьные оптимизации состояли в приведении подобных слагаемых и исключении слагаемых в (9) с нулевыми коэффициентами. Основная проблема, решавшаяся в рамках процедуры приведения подобных слагаемых — идентификация одинаковых (с точностью до порядка сомножителей) векторных частей в термах (9). Для решения этой проблемы было предложено использовать рекурсивную процедуру приведения векторного выражения к виду, минимальному в заданном лексикографическом порядке.
После проведения необходимых символьных вычислений и оптимизаций система генерирует исходных код на языке C++ для вычисления элементов локальных матриц и векторов в соответствии с полученными формулами. Помимо этого она строит выражения для нахождения интерполянтов функций J, E® и Hо, а также - для вычисления электрического и магнитного полей в тетраэдре.
2.4. Распределенные итерационные решатели
После декомпозиции задачи на подобласти и построения блоков матрицы на соответствующих узлах кластера, мы имеем распределенную расширенную СЛАУ
Ли = f.
Для решения таких СЛАУ широко используются методы Шварца. Одна итерация аддитивного метода Шварца в матричном виде может быть записана как (см. [10])
Ui+i = Ui + Е rTA-lRp(f - Ли),
p
где Rp — оператор сужения на подобласть Qp, p = 1,... ,P; здесь и далее подчеркивание для Ui мы опускаем. Матрицу Лр для подобласти p можно представить в виде Ap = RpARJ.
Проблема состоит в том, что для уравнения Гельмгольца метод Шварца может не сходиться. Однако можно заметить, что решение системы методом Шварца есть ни что иное, как решение методом простой итерации предобусловленной системы
В-1 Ли = В -1 f , В-1 = ^ RT A-1Rp. (10)
p
Такую предобусловленную систему можно решать итерационными методами в подпространствах Крылова. Однако метод Шварца имеет недостаток, заключающийся в том, что
количество итераций, необходимое для достижения заданной точности, растет с увеличением числа подобластей [10]. Для преодоления этого недостатка можно воспользоваться методом коррекции на грубой сетке [11] (Coarse Grid Correction, CGC).
Для этого введем операторы грубосеточного сужения Rc, определенный как оператор сужения на базис подпространства Wi, а также оператор продолжения Pc — R^. Тогда итерация метода аддитивного Шварца с грубосеточной коррекцией записывается как
ui+1/2 — ui + PcA- Rc( f Aui),
ui+1 = ui+1/2 + J2p Rp AplRp( f — A ui+1/2),
где A1 — RcAPc — матрица, полученная при использовании базисных функций порядка 1.
Мы будем решать предобусловленную систему вида (10), где предобуславливатель B-1 при наличии грубосеточной коррекции имеет вид
B-1 — Е RTA~lKp [I - APcA-1Rc] + PcA-1Rc. (11)
p
Поскольку метод коррекции на грубой сетке не очень чувствителен к точности аппроксимации A-1 [12], для приближенного решения систем вида A1V — q можно использовать итерационные методы в подпространствах следов, см. например [13]. Однако приближенное обращение матрицы A1 приводит к тому, что предобуславливатель B-1 становится переменным. В этом случае хорошим выбором для итерационного решения задач с таким пре-добуславливателем оказывается метод Flexible GMRES [10], допускающий использование переменных предобуславливателей.
Для решения задач в подобластях Apx — b можно использовать какой-нибудь прямой решатель для разреженных систем, например решатель PARDISO из библиотеки Intel MKL [6]. Проблема состоит в том, что время, требуемое PARDISO для разложения матриц, растет практически как O(N3/P3 ), поэтому такой метод подходит только для решения не очень больших систем на большом числе процессоров. Однако действие обратных матриц A-1 можно также вычислять приближенно, например методом FGMRES с предобуславли-вателем AMG, основанным на использовании иерархических базисов, предложенным в [13].
2.5. Методы визуализации электромагнитных полей
Визуализация электромагнитных полей является неотъемлемой частью численного эксперимента, поскольку позволяет проанализировать характерное распределение полей в расчетной области. В связи с этим, в рамках пакета Helmholtz3D был реализован модуль визуализации электромагнитных полей, основанный на технологии OpenGL.
Визуализация электромагнитных полей осуществляется следующим образом. Визуали-затор предоставляет возможность визуализировать одно из полей - электрическое либо магнитное, и позволяет переключаться между этими двумя режимами. Поля отрисовываются в каждом из тетраэдров в виде стрелки, сонаправленной с полем в данной точке, и имеющей длину, пропорциональную амплитуде поля в текущий момент времени. Для большей наглядности амплитуда поля кодируется цветом - с линейным переходом от желтого (максимально возможная амплитуда для текущего поля), до синего (близкая к нулю амплитуда). Отметим, что поскольку амплитуда поля представляется в виде двух векторов -действительной и мнимой компоненты, то в заданный момент времени t реальная ампли-
туда поля задается формулой
Р(I) = Рге СОё(ш1) + Ргт
где Рге и Ргт — векторы амплитуд действительной и мнимой части для соответствующего поля (электрического либо магнитного). Для примера на рис. 3 приведена визуализация электрического поля в волноводе с однородной средой. На рисунке хорошо видна структура электрического поля.
Рис. 3. Пример визуализации электрического поля в волноводе
3. Технологии параллельного программирования
_ __ о у'"' о
и оптимизации для систем с распределеннои и общей памятью
Параллельные распределенные алгоритмы в пакете Helmholtz3D построены на основе широко распространенного стандарта интерфейса обмена данными MPI (Message Passing Interface). При этом используемые в пакете итерационные решатели с предобуславливателя-ми типа аддитивного Шварца естественным образом ложатся на архитектуру современных вычислительных систем. Действительно, итерационное решение СЛАУ методом FGMRES требует 3 следующих типа операций.
К первому типу относятся векторно-векторные операции. Для имеющейся декомпозиции задачи на подобласти распараллеливание таких операций проводится естественным образом, за исключением операций вычисления скалярного произведения и нормы, для которых необходима точка синхронизации и обмен данными. Однако в интерфейсе MPI имеется требуемая функция MPI Allreduce, реализация которой оптимизируется поставщиком библиотеки MPI.
Для умножения матрицы на вектор решатель строит в каждой подобласти списки узлов, являющихся граничными для соседних подобластей. В дальнейшем решатель использует собранную информацию для пересылки частей вектора между процессами. Ключевым моментом в таком подходе является уменьшение объема пересылаемых данных — в соседние подобласти требуется пересылать только граничные коэффициенты векторов вместо того, чтобы пересылать вектор целиком. В итоге, общий объем пересылаемых данных в трехмерных задачах электромагнетизма для одного умножения матрицы на вектор будет пропорционален O(N2/3 log P) вместо O(N).
Применение предобуславливателей, основанных на методе Шварца, требует на каждой итерации внешнего алгоритма решения подзадач в подобластях, и, согласно (10), не требуют коммуникаций. Распараллеливание решателя в подобластях проводилось с использова-
нием средств OpenMP. Были реализованы параллельные процедуры умножения матрицы на вектор, а также параллельное решение разреженных треугольных систем для оператора сглаживания SSOR (Symmetric Successive Over-Relaxation, [10]) в предобуславливателе AMG. Распараллеливание треугольного решателя проводилось при помощи метода планирования вычислений по уровням: при решении системы вида Ux = у для каждой вершины вычисляется ее уровень li:
li = max {lj + 1 : Uiyj = 0, j > i} . (12)
После этого коэффициенты вектора x с одинаковым уровнем можно вычислять независимо.
Можно отметить ряд полезных оптимизаций решателя в подобластях. Для начала необходимо отметить его особенность: такое итерационное решение СЛАУ в подобластях оказывается bandwidth-limited, то есть существенно ограничено пропускной способностью шины данных системы, поскольку за одну итерацию алгоритма, фактически, требуется несколько раз прочитать всю матрицу системы (в том числе для решения двух треугольных систем), и несколько раз прочитать и записать различные векторы. Таким образом, реализованные оптимизации можно разделить на две категории: оптимизации доступа к памяти за счет переупорядочивания матрицы и оптимизации, направленные на повышение производительности алгоритма на NUMA-системах, заключающиеся в правильной инициализации областей памяти. Подробнее об этих оптимизациях см. [14].
Отметим, что помимо матрично-векторных операций, решатель в подобластях использует прямой метод PARDISO из Intel® MKL для самой «грубой» системы в предобуславливателе AMG. Дополнительное распараллеливание и оптимизация данной части решателя не требуется: PARDISO поддерживает OpenMP параллелизм как на этапе построения LU-разложения матрицы, так и на этапе поиска решения, а также использует высокоэффективные функции BLAS и LAPACK.
4. Численные эксперименты
Предлагаемые в работе алгоритмы были апробированы на серии методических и практических задач электромагнетизма.
Первой из задач является вычисление электромагнитного поля в параллелепипедаль-ном волноводе [0, а] х [0, b] х [0, с], а = 72 мм, b = 34 мм и с = 200 мм, см. рис. 4.
Рис. 4. Расчетная область для модельной задачи
Физические параметры среды: ^г = 1, ¿г = 1 — 0.И, ш = 6п ■ 109 Гц. При г = 200 мм поставлено краевое условие (3) с Е0 = еу $,т(жх/а), на остальной части границы — краевое
условие идеального проводника Eq = 0. Известен аналитический вид решения:
J._/ - \_У
sm y с
Сетка в данном случае была сгенерирована средствами пакета Helmholtz3D.
E = е
пх\ sin yz
(13)
Рис. 5. Расчетная область для задачи электромагнитного каротажа
Вторая из рассматриваемых задач — задача электромагнитного каротажа. Рассматриваемая расчетная область в виде куба с центром в начале координат и стороной 10 метров. Параметры цr = 1, ег = 1 во всей области. Проводимость среды а = 0.1 См/м, в среде имеется слой с а = 0.05 См/м. Координаты слоя по оси Oz: от -0.425 м до -0.275 м. В среде пробурена вертикальная цилиндрическая скважина радиуса 0.108 м с центром в начале координат в плоскости Oxy. В скважине имеется цилиндрическая каверна внешнего радиуса 0.118 м и положением по оси Oz от -0.0725 до 0.0725 м. Проводимость скважины и каверны а = 5 См/м. В скважину вставлен полый цилиндрический прибор радиуса 0.043 м, смещенный по оси x относительно центра скважины на -0.064, заполненный воздухом (а = 0). В приборе имеется одна генераторная и две приемные петли радиуса 0.0365 м при z = 0.5 м и при z = 0.0 и z = 0.1 м соответственно. По генераторной петле течет ток 1 А с частотой 14 МГц. Искомой является разность фаз ЭДС в приемниках.
Генерация неравномерной сетки проводилась при помощи утилиты NETGEN [7]. Количество тетраэдров сетки Nt = 490282, порядок СЛАУ без учета пересечений Nq = 8946864.
Тестирование проводились на двух системах: OpenMP решатель для систем с общей памятью (соответствует P = 1) — на системе Intel® Xeon X7560, 2.27 ГГц, 4 х 8 ядер, 64 ГБ памяти, гибридный решатель — на узлах HP БЬ2х220 G6 кластера НКС-30Т. Каждый узел кластера — Intel Xeon Е5540, 2.53 ГГц, 2 х 4 ядер, 16 ГБ памяти.
Для аппроксимации использовались базисные функции W3 третьего порядка.
В качестве критерия останова итерационного процесса было выбрано уменьшение нормы невязки в е раз: |г*| < е|/1. Для внутренних итераций FGMRES с предобуславливателем AMG ел = 0.2, для грубосеточной коррекции в подпространстве следов ес = 0.1. Решение задач в подобластях при обращении Ap проводилось с ер = 0.01. Внешние итерации проводились с точностью е = 10-7 за исключением сетки 116 х 55 х 320 для первой задачи, для которой использовалось е = 10-9 для получения Eh нужной точности.
В табл. 1-3 приведены результаты тестирования решателей на модельной задаче и задаче электромагнитного каротажа: Nq порядок системы без учета пересечений, SE — мак-
1 Используется более одного MPI процесса на узел (всего 64 узла).
Таблица 1
Масштабируемость гибридного кластерного решателя на модельной задаче
Сетка, N0 Количество подобластей (MPI процессов)
КЗ и 16 32 64 1281 2561 5121
58 х 28 х 160 Nb 705312 1168248 1777626 2599050 3928914 5804454
n 8 9 10 10 11 12
28519914 Пс 266 385 382 468 577 704
ПА 24 27 29 29 32 36
1 1 О 1 t 1.2e3 5.2e2 2.3e2 1.1e2 8.0e1 5.7e1
116 х 55 х 320 Nb 13035042 18153966
n 13 14
225335973 Пс N/A N/A N/A N/A 1102 1142
ПА 40 42
9,0 ■ 10-8 t 1.2e3 9.2e2
симальная относительная ошибка электрического поля, N — размерность пространства следов для Л\, п — количество внешних итераций распределенного решателя, пз и п — общее количество внешних и внутренних (в предобуславливателе ЛМО) итераций ОрепМР решателя, пс — общее количество итераций грубосеточной коррекции, па — максимальное количество внешних итераций метода РСМИЕВ в подобластях, Ь — время решения СЛАУ в секундах в экспоненциальном формате а.Ьеп = а.Ь • 10п. Результат решения задачи электромагнитного каротажа показал хорошее соответствие с результатом группы НГТУ, которая решала эту задачу другим методом — методом выделения поля источника, при этом искомая разность фаз ЭДС в приемниках совпала с ошибкой 1.8 %.
Таблица 2
Масштабируемость OpenMP решателя на задаче электромагнитного каротажа
Количество OpenMP потоков
1 2 4 8 16 32
П3 14 14 14 14 14 14
П2 33 33 33 33 33 33
t 7.5e3 4.1e3 2.2e3 1.3e3 8.8e2 6.1e2
Из таблиц видно, что предобуславливатели AMG и CGC успешно уменьшают количество внешних итераций до приемлемого значения. Количество итераций грубосеточной коррекции растет, но гораздо медленнее, чем линейно, за счет трехмерной декомпозиции расчетной области. Тем не менее, при дальнейшем росте размера СЛАУ может потребоваться геометрическое огрубление задачи для Л\, например [15], в результате чего мы получим многоуровневую схему грубосеточной коррекции.
1 Используется более одного MPI процесса на узел (всего 64 узла).
Таблица 3
Масштабируемость MPI решателя на задаче электромагнитного каротажа
Количество подобластей (MPI процессов)
4 8 16 32 64 1281 2561 5121
Nb 102918 236787 477672 803112 1346418 2038260 2873652 3972480
n 14 19 19 36 35 36 38 42
Пс 86 131 148 310 310 341 383 445
ПА 49 61 61 111 109 113 123 126
t 8.8e2 4.9e2 2.3e2 1.9e2 9.2e1 6.2e1 5.5e1 4.9e1
Заключение
В работе представлен пакет параллельных прикладных програм 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. Ingelstrom, P. A new set of H(curl)-conforming hierarchical basis functions for tetrahedral meshes / P. Ingelstrom // IEEE Transactions on Microwave Theory and Techniques. — 2006.
— 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. Schoberl, J. NETGEN — An advancing front 2D/3D-mesh generator based on abstract rules / J. Schoberl // Computing and Visualization in Science. — 1997. — Vol. 1, № 1. — P. 41-52.
8. Бутюгин, Д.С. Алгоритмы решения СЛАУ на системах с распределенной памятью в применении к задачам электромагнетизма / Д.С. Бутюгин // Вестник ЮУрГУ. Серия «Вычислительная математика и информатика». — 2012. — № 46(305). — С. 5-18.
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 deflation and coarse grid correction applied to porous media flow / R. Nabben, C. Vuik // SIAM Journal on Numerical Analysis. — 2004. — Vol. 42, № 4.
— P. 1631-1647.
13. Butyugin, D.S. Efficient 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 Scientific 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 fields 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 finite 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 field computation in the microwave devices and well logging problems.
Keywords: Maxwell’s equations, finite 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.
3. Ingelstrom P. A new set of H(curl)-conforming hierarchical basis functions for tetrahedral 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 finite 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. Schöberl J. NETGEN — An advancing front 2D/3D-mesh generator based on abstract rules. 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 deflation and coarse grid correction applied to porous media flow. SIAM Journal on Numerical Analysis. 2004. Vol. 42, No. 4. P. 1631-1647.
13. Butyugin D.S. Efficient 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 meAhods 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 Scientific Computing. 2005. Vol. 27, No. 5. P. 1669-1688.
Поступила в редакцию 15 февраля 2013 г.
