CC BY
105
АЛГОРИТМЫ РЕШЕНИЯ СЛАУ НА СИСТЕМАХ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ В ПРИМЕНЕНИИ К ЗАДАЧАМ ЭЛЕКТРОМАГНЕТИЗМА1
Д.С. Бутюгин
Рассматриваются различные аспекты моделирования гармонических электромагнитных полей на кластерах. Основная вычислительная сложность задачи заключается в решении систем линейных алгебраических уравнений (СЛАУ), возникающих в результате конечно-элементных аппроксимаций соответствующих краевых задач электромагнетизма элементами Неделека различных порядков. Рассмотрены эффективные и экономичные подходы к декомпозиции расчетной области и матрицы системы. Решение распределенных СЛАУ осуществляется итерационными методами в подпространствах Крылова с использованием аддитивного метода Шварца в качестве предобуславливателя. Для повышения эффективности алгоритмов итерации осуществляются в подпространствах следов. Реализованные решатели используют MPI для организации обмена данными. Решение систем в подобластях осуществляется при помощи прямого решателя PARDISO из библиотеки Intel® MKL. Результаты серии численных экспериментов на модельных и практических задачах демонстрируют эффективность предлагаемых алгоритмов.
Ключевые слова: уравнения Максвелла, итерационные алгоритмы, методы декомпозиции подобластей, аддитивный метод Шварца.
Введение
Многие актуальные задачи вычислительной математики требуют решения разреженных СЛАУ высоких порядков. Задача вычисления трехмерного электромагнитного поля в частотной области возникает при расчетах различных волновых устройств, решении задач геоэлектроразведки, таких как электромагнитного каротажа, и других. Требования к точности получаемого решения задач с разномасштабными объектами приводят к необходимости построения сеток с числом конечных элементов до 105 ^ 107 и порядком получаемых после аппроксимации СЛАУ до 107 ^ 109. Решение таких СЛАУ невозможно без использования вычислительной мощности кластеров. Это требует подходящих алгоритмов решения систем алгебраических уравнений для машин с распределенной памятью.
В данной работе рассматриваются алгоритмы решения СЛАУ, основанные на аддитивном методе Шварца. Рассмотрено два подхода к декомпозиции задачи. Первый из подходов основан на геометрическом разбиении сетки на связанные подобласти и позволяет построить декомпозицию с пересечениями подобластей. Второй подход основан на переупорядочивании графа матрицы системы. Этот способ представляет из себя вариант одно-направленного разбиения графа и приводит к СЛАУ блочнотрехдиагонального вида. Отличительной особенностью предлагаемых методов является низкая асимптотическая вычислительная сложность алгоритмов, в отличие от вариантов метода вложенных сечений (например, реализованного в библиотеке METIS [1]).
1 Статья рекомендована к публикации программным комитетом международной научной конфе-
ренции «Параллельные вычислительные технологии 2012»
Решение полученных систем уравнений осуществляется методами в подпространствах Крылова [2], такими как обобщенный метод минимальных невязок (GMR.ES). В качестве предобуславливателя выступает аддитивный метод Шварца. Использование такого подхода, в отличие от использования метода Шварца напрямую, позволяет решить более широкий класс задач и повысить скорость сходимости итерационного процесса.
Структура данной работы следующая. В разделе 1 описывается математическая постановка задачи, в разделе 2 — дискретная постановка. В разделе 3 описаны предлагаемые методы декомпозиции. Раздел 4 содержит описание используемых итерационных алгоритмов. В разделе 5 приведены результаты численных экспериментов. Наконец, в последнем разделе обсуждаются полученные в работе результаты.
1. Математическая постановка задачи
Электромагнитное поле с гармонической зависимостью от времени в случае линейности задачи при отсутствии внешнего заряда и магнитной проводимости может быть описано следующей формой системы уравнений Максвелла:
Ух Е = —гш^Н, V ■ (ег Е) = р/е 0, (1)
Ух Й = гшеЕ + 3, V ■ (^г Н) = 0,
где электрическая и магнитная проницаемости имеют вид
е = еоег — 1(Ге/ш, ^ = ^о^г — %от/ш.
Здесь г — мнимая единица, Е и Н — векторы напряженности электрического и магнитного полей соответственно, ш — круговая частота решения, . и р — плотности внешнего тока и объемного заряда соответственно, е0 и ^0 — диэлектрическая и магнитная проницаемости вакуума, ег и ^г — относительная диэлектрическая и магнитная проницаемости среды, а ае и ат — электрическая и магнитная проводимости.
Используя предположения о линейности задачи (то есть независимости параметров среды от электромагнитного поля), а также отсутствие объемного внешнего заряда р и магнитной проводимости ат = 0, уравнения (1) можно привести к форме комплексного векторного уравнения Гельмгольца
У х V х Е^ — к^єгЕ — —.
Здесь к0 = шл/е0^0, Z0 = \/^0/е0 и ег = е/е0. Параметры среды полагаются кусочно-постоянными. Мы будем предполагать, что к0 не является Максвелловским собственным числом, то есть рабочая частота ш не является резонансной. Стоит отметить, что в случае, если хотя бы в одной подобласти ае > 0, то рабочая частота может быть любой, поскольку в системе отсутствуют резонансы [3].
Решение ищется в области П с границей 5 = 51 и Б2, на каждой из частей которой поставлено одно из следующих граничных условий:
п х Е
— п х Е0, п х Н
Sl
— 0, (3)
S2
где п — внешняя нормаль к границе. При Е0 = 0 первое из условий соответствует идеальному электрическому проводнику, при ЕЙ0 = 0 — волновому входу, а второе
условие соответствует идеальному магнитному проводнику. Полагая область П состоящей из подобластей П = У П, в каждой из которых физические параметры среды являются константными, введем на каждой из внутренних границ Г с нормалью п условия сопряжения
п ■ (ёгЕг - ё^) = 0, п х (Ег - Е2) = 0,
п ■ (ргЙг - ^2^2) = 0, п х (Я 1 - Н2) = 0.
Вводятся стандартные соболевские пространства и подмножества:
Йг°* = {ф е [Ь2(П)]3 : Ух ф е [£2(П)]3} ,
Й0Г<Л = { ф е Йго1 : п х ф |51 = 0} , Щ°* = { ф е Йг<* : п х Я|31 = п х £0} .
Вариационная формулировка уравнения (2) с краевыми условиями (3) в форме Галёркина имеет следующий вид (см. [4]): найти такое Е е Й3°ь, что для всех ф е Й0° выполнено
! ^У х Е^ ■ (у х ф^ dП - к^ J £г (^Е ■ ф^ dП = —ik0Z0 ^ 3 ■ ф(5)
2. Дискретная постановка задачи
В работе используются неструктурированные тетраэдральные сетки. Рассматривается соответствующее разбиение расчетной области П на непересекающиеся тетраэдральные элементы П = и Пт. В каждом из тетраэдров вводятся базисные функции, соответствующие его степеням свободы. Пусть Щг — конечномерное пространство базисных функций порядка не выше I, конформное Йго1. В работе рассматриваются иерархические базисные функции, предложенные в [5]:
Щ = Щ ФУУ2 0 ... ФУЩ,
Щ = Аг, Щ = А ФУУ*,
где УЩ — инкрементальное подпространство с базисными функциями порядка i, V — инкрементальное подпространство со скалярными базисными функциями, конформными Йг, а А — инкрементальное подпространство с роторными базисными функциями.
Подпространства VI,0 и Нг,0 вводятся естественным образом, а подмножество
— как множество функций, у которых коэффициенты разложения п3 по базисным функциям с ненулевым следом на Бг принимают фиксированные значения, соответствующие краевому условию (3).
Приближенное решение Ен будем искать в виде
Ен = £ + £ „#». (6)
3 г
Вводятся матрицы соответствующих билинейных форм:
Мг,з = I £г (ф° ■ ф°) dП, Бг,з = I (у х ф°) ■ (у х ф°) dП,
где базисные функции щ ,'ф0 е Нг,0. Вектор правой части определяется как
/г = ° Н ■ ЩМП + ^ п3 ! (р-1 (у х 3 ■ (у х - к1егН3 ■ dП.
' ° 3 °
Тогда итоговая система принимает вид
[5 - к2вЫ] п = /. (7)
Для вычисления элементов матрицы и вектора правой части можно воспользоваться поэлементной технологией сборки [6], заменив интегрирование по расчетной области П суммой интегралов по каждому из тетраэдров и вычислением в каждом из тетраэдров локальных матриц и векторов.
3. Методы декомпозиции на подобласти
Для решения задачи на системах с распределенной памятью требуется распределить данные между узлами системы, по возможности минимизируя объем коммуникационных данных. Основной проблемой здесь является то, что алгоритмы, генерирующие разбиения высокого качества, сами требуют больших вычислительных ресурсов. Поэтому, как правило, приходится идти на компромисс и выбирать алгоритмы, позволяющие получить достаточно хорошее разбиение за приемлемое время.
3.1. Геометрическая декомпозиция расчетной области
Предлагаемый в работе алгоритм геометрической декомпозиции расчетной области является упрощенным вариантом алгоритма построения BSP-дерева (binary space partitioning, двоичное разбиение пространства). Метод двоичного разбиения — это метод рекурсивного разбиения пространства на выпуклые множества гиперплоскостями [7]. Получающаяся в результате структура данных находит широкое применение в различных областях компьютерной графики для проверки пересечений, определения прямой видимости объектов и др. Стоит отметить, что построение оптимального BSP дерева представляет из себя очень трудоемкую задачу, поэтому, в большинстве случаев, ограничиваются субоптимальными деревьями.
В рамках данной работы предлагается использовать следующее упрощенное двоичное дерево. Упрощение состоит в том, что секущие плоскости предлагается проводить ортогонально осям координат. При этом мы будем каждый раз разделять имеющееся множество тетраэдров на две приблизительно равные части. Алгоритм построения разбиения будет основан на методе обработки событий. Подготовительная часть алгоритма такая:
• для х, у и z: спроектировать тетраэдры на соответствующую ось; для каждого тетраэдра получаем отрезки, соответствующие его началу и концу;
• для х, у и z: сгенерировать отсортированные по координате события — начало и конец каждого из тетраэдров.
В результате мы получим 3 отсортированных списка событий, для осей Ox, Oy и Oz соответственно.
Далее к этим спискам мы рекурсивно будем применять следующую процедуру. Будем последовательно обрабатывать события для каждой из осей, что соответствует движению секущей плоскости в сторону увеличения соответствующей координаты.
При этом мы будем подсчитывать число тетраэдров в каждом из полупространств. Тетраэдры, пересекающиеся и касающиеся плоскости сечения, будем относить к обеим частям одновременно. Схема алгоритма следующая:
• Для x, у и z:
— left_count ^ 0, right_count ^ T, где T — число тетраэдров.
— Для каждой координаты последовательно: обработать события
* Обработать начала тетраэдров: для каждого left_count ^ left_count +
1.
* Проверить left_count и right_count на оптимальность.
* Обработать концы тетраэдров: right_count ^ right_count — 1.
— Запомнить оптимальное разбиение для оси.
• Выбрать наилучшую ось для разбиения.
• Построить разбиение: для x, у и z:
— Линейным проходом разбить список событий для соответствующей оси на 2 в соответствии с разбиением. При этом упорядоченность сохранится.
• Рекурсивно запустить алгоритм для двух полученных частей.
Алгоритм останавливает рекурсивное разбиение тогда, когда количество тетраэдров в обрабатываемой части становится достаточно мало, либо когда в результате деление количество тетраэдров в одной из частей оказался равен исходному их числу.
Проверка на оптимальность для каждой из осей может быть следующей: выбирается положение плоскости, при котором размеры левой и правой частей максимально близки друг к другу. Наилучшая из осей выбирается по минимуму min(left_count, right_count) для оптимального положения плоскости, что позволяет выбрать минимальный из разрезов подпространства и уменьшить коммуникации между подобластями. Отметим, что «толщину» слоя пересечения легко варьировать. Для этого достаточно осуществить проход по событиям двумя указателями вместо одного. Один из указателей при этом будет соответствовать началу пересечения подобластей, второй — концу. При этом легко поддерживать как нужную геометрическую толщину пересечения подобластей, так и минимальное число тетраэдров в пересечении.
Легко видеть, что асимптотическая сложность алгоритма равна O(TlogT), где T — исходное число тетраэдров, при условии, что каждый раз в пересечении оказывается небольшое число тетраэдров и части делятся примерно поровну. Действительно, сложность сортировки событий в начале алгоритма равна O(T log T). Далее, если число тетраэдров растет несущественно по сравнению с величиной T (например, конечное число тетраэдров есть O(T)), то на каждом уровне дерева требуется линейный проход по всем событиям для каждого из узлов дерева этого уровня, суммарное количество которых есть O(T). При этом число подобластей при переходе на уровень ниже увеличивается в 2 раза, поэтому общее число уровней не превзойдет O(log T). Отсюда получаем заявленную сложность.
Замечание 1. Алгоритм осуществляет разбиение сетки на подобласти с пересечениями, причем никакие две подобласти не имеют протяженной общей границы кроме, возможно, части внешней границы расчетной области.
Доказательство этого факта следует из следующего тривиального наблюдения: два тетраэдра с общей гранью невозможно разбить плоскостью так, чтобы один тет-
раэдр лежал по одну сторону от нее, второй по другую и плоскость не касалась хотя бы одного из тетраэдров. При этом, как следует из алгоритма, тетраэдры, касающиеся либо пересекающиеся с плоскостью сечения, помещаются в пересечение подобластей. Данный факт оказывается очень важным для дальнейшего решения задачи методом декомпозиции на подобласти.
Рассмотрим теперь постановку задачи нахождения поля для подобластей. Пусть имеется разбиение расчетной области П на подобласти П = У Пр. Будем искать решение задачи (5) на области П путем решения задачи (5) в каждой из подобластей Пр. На внутренней границе Г каждой из подобластей Пр будем ставить условия Дирихле из (3). Эти условия будут определяться из значений поля на внутри подобластей, соседних с Пр, и, на конечно-элементном уровне, из коэффициентов разложения поля Е по базисным функциям, соответствующим граничным ребрам и граням. То есть, требуется найти поля в каждой из подобластей так, чтобы они были согласованы друг с другом по краевым условиям.
Теорема 1. Решение задачи для подобластей для предлагаемого алгоритма декомпозиции существует. Оно единственно, если частота излучения не совпадает с резонансными частотами полостей, образованных подобластями и их пересечениями.
Доказательство. Ясно, что решение исходной задачи на всей области П удовлетворяет задаче на подобластях. Единственность следует из следующей леммы.
Лемма 1. Пусть поле Е является решением задачи для подобластей. Тогда для любых Пр и Пр/, имеющих непустое пересечение, значения Е в общих точках Пр и Пр/ совпадают.
Доказательство напрямую следует из того, что для пересечения подобластей мы имеем краевую задачу, которая, если частота не резонансная для данной задачи, имеет единственное решение. Поэтому решение задачи для подобластей однозначно определяет электрическое поле в любой точке исходной области П.
Доказательство теоремы 1 завершается замечанием того, что
• для любой внутренней точки области П существует замкнутая окрестность, что эта точка вместе с окрестностью целиком лежит внутри какой-либо подобласти Пр;
• для любой базисной функции из Ніоо ее носитель целиком лежит внутри какой-нибудь подобласти Пр;
Отсюда следует, что электрическое поле, являющееся решением задачи для подобластей, удовлетворяет и вариационной формулировки (5) в непрерывном случае, и дискретному аналогу этой формулировке в конечно-элементном случае. Однако решение вариационной формулировки единственно при указанных ограничениях на параметры среды и частоту излучения [3], поэтому единственным оказывается и решение задачи для подобластей. □
3.2. Алгебраическая декомпозиция
Другим вариантом декомпозиции задачи на подобласти, рассматриваемым в рамках работы, является модифицированный алгоритм одно-направленного разбиения
графа матрицы системы [2]. Алгоритм строит разбиение неизвестных системы на множества — алгебраические подобласти П^. В дальнейшем индекс A в мы будем опускать.
Первым шагом алгоритма является нахождение псевдо-периферийной вершины v. Поиск такой вершины осуществляется следующим образом [2].
1. Выбирается произвольная вершина v. Псевдо-диаметр графа D полагается равным 0.
2. Запускается обход в ширину по графу, начиная от вершины v.
3. Находится любая максимально удаленная вершина V.
4. Если расстояние d(v,v') больше D, то v ^ V, D ^ d(v,v'), переход на шаг 2.
Информация с последнего обхода графа в ширину используется для дальнейшего
построения разбиения. Для начала все вершины разбиваются во «фронты» Fk — множества вершин, равноудаленных от v на расстояние к. Далее фронты объединяются в алгебраические подобласти по следующему принципу.
С помощью бинарного поиска мы будем максимизировать размер минимальной подобласти. Проверка того, что при данном минимальном размере подобласти можно получить необходимое число подобластей, может быть осуществлена жадным алгоритмом [8]. Пример кода, определяющего максимальное число подобластей, представлен на рис. 1.
int get_max_domains(int k, int front_size[], int min_size) { int domains = 0, domain_size =0, i; for(i =0; i < k; i++) {
domain_size += front_size[i]; // Размер в вершинах if(domain_size >= min_size) {
++domains; domain_size = 0;
}
}
// Оставшиеся вершины отходят к последней подобласти return domains;
}
Рис. 1. Алгоритм определения максимального числа подобластей
Представленный алгоритм генерирует подобласти без пересечений. Однако получить разбиение с пересечениями довольно легко — достаточно при начале генерации следующей подобласти отступить от границы предыдущей внутрь необходимое число фронтов. Графическое представление объединения фронтов в подобласти с пересечениями изображено на рис. 2, где фронты изображены вертикальными линиями.
\
/
Рис. 2. Объединение фронтов матрицы в подобласти
Or
O
p+1
Обозначив 1р и гр — границы подобласти Пр, мы имеем разбиение на подобласти
Гр
^р = 1^) Рк, к—1р
Необходимо отметить, что ребра в графе матрицы могут соединять вершины либо одного, либо двух соседних фронтов. Это означает, что каждая подобласть Пр может граничить одновременно максимум с двумя подобластями — Пр_1 и Пр+1 (при наличии пересечений всегда можно считать, что два граничных с Пр фронта Flp_1 и +1 принадлежат соседним подобластям Пр-1 и Пр+1 соответственно). Тогда последовательно занумеровав переменные в каждой из подобластей, а в каждой из подобластей — последовательно в каждом из фронтов, мы получим следующую блочнотрехдиагональную СЛАУ:
■ Рі и и1 ' /і '
Ьі п2 /2
••• ир-1 1 ./р 1
Ь р і р 1 ир
При этом если разбиение исходной СЛАУ выполнялось с пересечениями, то в итоговой СЛАУ окажется по несколько переменных, соответствующих вершинам в пересечении.
4. Алгоритмы решения СЛАУ
После декомпозиции задачи на подобласти и построения блоков матрицы на соответствующих узлах кластера, мы имеем распределенную СЛАУ
Аи = /.
Для решения таких СЛАУ широко используются методы Шварца. Одна итерация аддитивного метода Шварца в матричном виде может быть записана как (см. [2])
иг+1 = иг + ^ ] ДрА_ 1^р(/ — Аиг), р
где Яр — оператор сужения на подобласть Пр. При этом, на самом деле, итерировать достаточно только переменные, соответствующие внутренним границам подобластей, поскольку зная граничные значения, легко определить значения внутри подобластей, решив соответствующую задачу в подобласти. То есть
ир = А-1Др(/ — АЯТ иГ), (8)
где иг — вектор граничных переменных, а Д — оператор сужения на граничные переменные. Подпространство граничных переменных принято называть подпространством следов. Тогда итерация метода Шварца в подпространстве следов переписывается как
иГ+1 = иГ + ДТА_1 Др(/ — А^ТиГ)-
р
Эту итерацию можно представить в виде иг+1 = Б (иг) = Тиг + д, где
Т = I — Я ^ ДТА-1ЯрАЯт, д = Я ^ Ятр А-1 Яр/. рр Сходимость данного метода имеется при спектральном радиусе р(Т) < 1. В то время как для эллиптических уравнений это условие выполнено, для уравнения Гельмгольца это не так. Однако в этом случае можно рассмотреть решение системы
[I — Т] иг = д. (9)
Теорема 2. Система (9) является совместной. В случае геометрической декомпозиции при выполнении условий теоремы 1 оно единственно.
Доказательство. Действительно, система (9) эквивалентна системе
ЯУ^ Я'ТА_1ЯрАи = Я^2 ЯТА-1 Яр/, рр поэтому решение иГ = Яи, где и — решение Аи = / удовлетворяет ей. С другой стороны любое решение (9) является неподвижной точкой итерации Шварца. В случае геометрической декомпозиции на подобласти у неподвижной точки итерации Шварца есть простая интерпретация. Это вектор, задающий поле в подобластях, согласованное по граничным условиям. Однако при выполнении условий теоремы 1 последняя утверждает, что такое поле единственно. Значит, единственна и неподвижная точка итерации Шварца. □
Можно отметить, что решение системы методом Шварца есть решение предо-бусловленной системы М-1Аи = М-1 / с М-1, связанным с А и Т соотношением Т = I — М-1 А, откуда М-1 = [I — Т]А-1 и предобусловленная система переходит в систему (9). Эту систему можно решать и итерационными методами в подпространствах Крылова, например методом обобщенных минимальных невязок (GMR.ES) [2]. Для методов в подпространствах Крылова не требуется явный вид матрицы системы, а достаточно только операции умножения матрицы на вектор. Запишем результат умножения I — Т на вектор V через результат одной итерации метода аддитивного Шварца Б(V):
[I — Т] V = V — ^ = V — Б^)+ д, д = Б(0).
Таким образом, для решения СЛАУ (9) методом GMRES достаточно реализовать обычный метод аддитивного Шварца. Для решения задач в подобластях Арх = Ь можно использовать какой-нибудь прямой решатель для разреженных систем. В работе использовался решатель PARDISO из библиотеки Ш;е1 МКЬ [9]. Такой подход оказывается достаточно эффективным, поскольку на каждой итерации требуется решать системы с одними и теми же матрицами, и Ш-разложение матриц Ар можно посчитать только один раз. Кроме того, такой подход оказывается существенно эффективнее простого использования решателя PARDISO для исходной системы, поскольку время, требуемое PARDISO для разложения матрицы, растет практически как О(М3), где N — порядок системы. Дополнительно эффективность достигается
за счет использования подпространства следов, так как в этом случае существенно уменьшается объем коммуникационных затрат, а также объем памяти, требуемый для хранения базиса подпространства Крылова в методе GMRES, поскольку хранить и пересылать требуется только граничные переменные. После того, как решение системы найдено, искомое поле Е восстанавливается в каждой из подобластей.
5. Численные эксперименты
5.1. Модельная задача
Рис. 3. Расчетная область для модельной задачи
В качестве одного из тестов рассматривалась модельная задача с расчетной областью в виде волновода с линейными размерами a = 72, b = 34, c = 200 мм (см. рис. 3. Параметры задачи: рr = 1, ег = 1 — 0.1 i, частота ш = 6п • 109 Гц. Граничные условия: на грани z = 200 задавалось условие
Е0 х n = ey sin(n x/a),
на остальных гранях — условие металлической стенки (Е0 = 0). Аналитическое решение
^ . /пх\ sinyz
Е = ey sin — ---------.
V a / sin yc
При решении СЛАУ критерием остановки итерационного процесса служило выполнение условия (r, r) < e2(f, f), где f — вектор правой части системы, r = f — Au — вектор невязки. В проведенных экспериментах е было выбрано равным 10-7.
Для построения сетки расчетная область делилась по каждому из измерений на некоторое число шагов, и каждый из получающихся параллелепипедов делился на 6 тетраэдров. В тестах использовалась квазиравномерная сетка с разбиениями по каждой из координат с некоторым шагом h. Были рассмотрены сетки с разбиениями 4 х 2 х 10, 8 х 4 х 20 и 15 х 7 х 40. Для аппроксимации использовались базисные функции второго порядка.
Тесты проводились на двух системах: системе Intel Core2 Duo E6550 @ 2.33 ГГц, 2 ядра, с оперативной памятью объемом 4 ГБ, и на узлах HP БЬ2х220 G6 кластера НКС-30Т. Каждый узел — Intel Xeon Е5540, 2.53 ГГц, 2x4 ядер, 16 ГБ памяти.
В табл. 1 приведены результаты тестирования решателя при использовании геометрической декомпозиции. В таблице nproc — число MPI процессов (и подобластей), N — порядок системы для nproc=1, Nb — размерность пространства следов, Np — максимальный порядок СЛАУ для подобластей, n — число итераций решателя GMRES,
5E — относительная ошибка для поля E в барицентрах тетраэдров. Тесты для сетки 4 х 2 х 10 при 16, 32 и 64 подобластях не проводились, так как при этом не хватает тетраэдров сетки для построения нужного числа подобластей. Тесты проводились на двухъядерном Core2 Duo, поэтому время работы решателей не приводится.
Таблица 1
Решение модельной задачи с геометрической декомпозицией
Сетка N 8E nproc
2 4 8 16 32 64
4 x 2 x l0 2392 3,1e-02 n 8 11 22 — — —
Nb 136 408 2656 — — —
Np 1408 916 670 — — —
8 x 4 x 20 21664 7,3e-03 n 9 14 23 40 53 70
Nb 592 1776 4556 15764 40396 105200
Np 11782 6292 4522 2610 1834 1058
l5 x Т x 40 149874 2,1e-03 n 10 17 27 50 64 80
Nb 2012 6036 14084 34560 67872 143256
Np 78206 40486 21626 12818 8346 5478
В табл. 2 приведены результаты тестирования для случая алгебраической декомпозиции. Параметр пересечения подобластей был установлен в 4, то есть подобласти пересекались, если это было возможно, по 4 фронтам. Из таблицы видно, что при декомпозиции матрицы на слишком большое число подобластей число итераций начинает быстро расти. В этих случаях на каждую подобласть приходится всего по нескольку фронтов.
Таблица 2
Решение модельной задачи с алгебраической декомпозицией
Сетка N 8E nproc
2 4 8 16 32
4 x 2 x l0 2392 3,1e-02 n 6 12 49 — —
Nb 378 980 2261 — —
Np 1542 1334 796 — —
8 x 4 x 20 21664 7,3e-03 n 9 13 21 258 —
Nb 1540 4678 10708 23052 —
Np 12420 9232 6306 3200 —
l5 x Т x 40 149874 2,1e-03 n 13 18 25 38 517
Nb 5210 16646 37670 78446 162720
Np 81132 49618 31168 22678 11282
5.2. Задача электромагнитного каротажа
Вторая из рассматриваемых задач — задача электромагнитного каротажа. Расчетная область представляет из себя куб с центром в начале координат со стороной 10 метров. Куб заполнен средой с а = 0,1 См/м, в среде имеется слой с а = 0, 05
См/м. Координаты слоя по оси Oz: от —0, 425 м до —0, 2Т5 м. В среде пробурена вертикальная цилиндрическая скважина радиуса 0, l08 м с центром в начале координат в плоскости Oxy. В скважине имеется цилиндрическая каверна внешнего радиуса
0, ll8 ми положением по оси Oz от —0, 0Т25 до 0, 0Т25 м. Проводимость скважины и каверны о = 5 См/м. В скважину вставлен полый цилиндрический прибор радиуса
0,043 м, смещенный по оси x относительно центра скважины на —0, 064. Проводимость прибора равна нулю. В приборе имеется одна генераторная петля при z = 0,5 м и две приемные петли при z = 0, 0 и z = 0, l м. Радиус петель — 0, 0365 м. По генераторной петле течет ток 1 А с частотой 14 МГц. Искомой является разность фаз ЭДС в приемниках.
Генерация неравномерной сетки проводилась при помощи утилиты NETGEN v4.9.13 [10]. В результате была получена сетка с 71892 узлами и 388836 тетраэдрами.
Результаты решения задачи на кластере с использованием алгебраической декомпозиции представлены в табл. 3. В ней дополнительно указаны времена tfact — максимальное время факторизации матриц Dp в подобластях и ttot — общее время работы решателя. Результаты тестирования для nproc < 4 не приведены, так как в этом случае PARDISO не хватает 16 Гб памяти на узле для хранения LU факторов матрицы, а режим Out-of-Core с хранением факторов на диске является слишком медленным.
Таблица 3
Решение задачи каротажа с алгебраической декомпозицией матрицы
nproc N Nb tfact, с ttot, с n
4 2398750 371126 8,52e+02 9,49e+02 21
8 2398750 833744 3,30e+02 4,26e+02 28
Результат решения задачи показал хорошее соответствие с результатом группы НГТУ, которая решала эту задачу другим методом — методом выделения поля источника, при этом искомая разность фаз ЭДС в приемниках совпала с ошибкой 1.8%. Стоит отметить также, что для данной задачи обычные итерационные алгоритмы в подпространствах Крылова, примененные к СЛАУ для всей подобласти, показывают отсутствие сходимости.
Заключение
В работе рассматривается задача моделирования электромагнитного поля в частотной области. Предложено два эффективных алгоритма геометрической и алгебраической декомпозиции задачи, причем последний может применяться и для решения других задач. Предлагаемые итерационные алгоритмы в подпространствах следов с предобуславливателем Шварца требуют невысоких коммуникационных затрат и хорошо подходят для вычислительных систем с распределенной памятью. Реализация алгоритмов на языке С+—Н с использованием средств МР1 продемонстрировала хорошую производительность на методических и практических задачах. Предложенные решатели позволили решить задачу электромагнитного каротажа с высокой точностью.
Работа выполнена при поддержке РФФИ (грант 11-01-00205).
Литература
1. Karypis, G. A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular Graphs / G. Karypis, V. Kumar // SIAM Journal on Scientific Computing. - 1999. -Vol. 20, № 1. - P. 359-392.
2. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition. / Y. Saad -Society for Industrial and Applied Mathematics, 2003.
3. Monk, P. Finite Element Methods for Maxwell’s Equations. / P. Monk - Oxford University Press, 2003.
4. Соловейчик, Ю.Г. Метод конечных элементов для решения скалярных и векторных задач / Ю.Г. Соловейчик, М.Э. Рояк, М.Г. Персова - Новосибирск: Изд-во НГТУ, 2007.
5. 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.
6. Ильин, В.П. Методы и технологии конечных элементов. / В.П. Ильин - Новосибирск: Изд-во ИВМиМГ СО РАН, 2007.
7. 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.
8. Кормен, Т. Алгоритмы: построение и анализ / Т. Кормен, Ч. Лейзерсон, Р. Ривест
- М., МЦНМО: БИНОМ. Лаборатория знаний, 2004.
9. Intel. The Flagship High-Performance Computing Math Library for
Windows*, Linux*, and Mac OS* X. Intel® Math Kernel Library from Intel. URL: http://software.intel.com/en-us/articles/intel-mkl/ (дата обращения:
22.01.2012)
10. 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.
Дмитрий Сергеевич Бутюгин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, аспирант, Новосибирский государственный университет, dm.butyugin@gmail.com.
ALGORITHMS OF SLAES SOLUTION FOR THE SYSTEMS WITH DISTRIBUTED MEMORY APPLIED TO THE PROBLEMS OF ELECTROMAGNETISM
D.S. Butyugin, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation)
Paper presents various aspects of harmonic electromagnetic fields simulation on clusters. The major computational complexity comes from the solution of the systems of linear algebraic equations (SLAEs) arising from the approximations of corresponding electromagnetic boundary value problems by Nedelec elements of various orders. Effective and efficient approaches to the decomposition of the computational domain and the matrix of the system are considered. Distributed SLAEs are solved using iterative Krylov subspace methods preconditioned by additive Schwarz method. In order to increase the effectiveness of the algorithms iterations are performed in the trace space. Implementation of the solvers is based on MPI for data transfers. The solution of the systems in subdomains is performed by PARDISO direct solver from Intel® MKL library. Numerical experiments results on a series of model and real-life problems show the effectiveness of the presented algorithms.
Keywords: Maxwell equations, iterative algorithms, domain decomposition methods, additive Schwarz method.
References
1. Karypis G., Kumar V. A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular Graphs. SIAM Journal on Scientific Computing. 1999. Vol. 20, № 1. P. 359392.
2. Saad Y. Iterative Methods for Sparse Linear Systems, Second Edition. Society for Industrial and Applied Mathematics, 2003.
3. Monk P. Finite Element Methods for Maxwell’s Equations. Oxford University Press, 2003.
4. 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.
5. 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, № 1. P. 160-114.
6. Il’in V.P. Metody i tekhnologii konechnykh elementov [Methods and Technologies of Finite Elements]. Novosibirsk, ICM&MG SBRAS Publ., 2007.
7. Fuchs H., Kedem Z.M., Naylor B.F. On Visible Surface Generation by A Priori Tree Structures ACM Computer Graphics. 1980. Vol. 14, № 3. P. 124-133.
8. Cormen T., Leiserson C., Rivest R. Introduction to Algorithms. MIT Press, 2001.
9. Intel (R) Math Kernel Library from Intel:
URL: http://software.intel.com/en-us/articles/intel-mkl/
10. Schoberl J. NETGEN — an Advancing Front 2D/3D-mesh Generator Based on Abstract Rules Computing and Visualization in Science. 1997. Vol. 1, № 1. P. 4152.
Поступила в редакцию 14 марта 2012 г.
