CC BY
39
Вычислительные технологии
Том 18, № 2, 2013
Параллельные конечноэлементные вычислительные схемы в задачах геоэлектрики
М.И. Эпов1, Э.П. Шурина1'2, Д. А. Архипов1'2 1 Институт нефтегазовой геологии и геофизики СО РАН, Новосибирск, Россия 2Новосибирский государственный технический университет, Россия e-mail: shurina@online.sinor.ru, d_arhipov@list.ru
Решение задач геоэлектрики в сложнопостроенных средах с контрастными электрофизическими свойствами приводит к необходимости разработки стабильных, реализуемых на современных параллельных суперкомпьютерах вычислительных схем. В работе приводятся результаты моделирования трёхмерных электромагнитных полей методом декомпозиции области и мультипликативным методом Шварца, исследуется влияние на результаты вычислений меры налегания подобластей.
Ключевые слова: векторный метод конечных элементов, декомпозиция области, мультипликативный метод Шварца.
Введение
Задачи геоэлектрики можно охарактеризовать следующим образом: область моделирования — геометрически сложная структура с внутренними границами, разномасштабными фрагментами и контрастными электрофизическими свойствами этих фрагментов; временные или частотные характеристики систем источник — приемник, их расположение и геометрические особенности.
Необходимость разработки эффективных вычислительных схем для интерпретации результатов измерений, полученных в наземной или скважинной электроразведке, приводит к возрастанию требований к процедурам многовариантного прямого моделирования, т. е. к применению эффективных алгоритмов решения систем уравнений Максвелла в трёхмерных сложнопостроенных областях. В течение последних десяти лет для таких расчётов используется векторный метод конечных элементов (ВМКЭ) [1, 2]. Однако особенности гс^-оператора (наличие нуль-ядра [3]) и ассоциированная с ним знаконеопределённость матрицы конечноэлементной системы линейных алгебраических уравнений (СЛАУ), характерная при гармонической зависимости от времени (частотная область) на частотах от десятков килогерц до десятков мегагерц в средах с разрывными электрофизическими характеристиками, сложными межфрагментарными границами, разномасштабными фрагментами, обусловливают создание адаптированных к этому классу задач алгоритмов, ориентированных на применение современных суперкомпьютеров. Одним из таких направлений являются методы на основе декомпозиции области — геометрический или алгебраический варианты — и специальное предо-бусловливание таких СЛАУ, учитывающее особенности векторной конечноэлементной аппроксимации на адаптивных тетраэдральных сетках. В настоящей работе рассматривается параллельный алгоритм моделирования трёхмерного электромагнитного поля,
базирующийся на декомпозиции области с налеганием, исследуются варианты декомпозиции области и анализируется влияние налегания подобластей на вычислительную эффективность алгоритма.
1. Математическая модель
Система уравнений Максвелла при гармонической зависимости решения от времени (вш1) имеет вид [3]
го1Е = —гшБ, (1)
гоШ = гшБ + аЕ + (2)
а1у Б = 0, (3)
а1у Б = р, (4)
Б = еЕ, (5)
где Б = ^Ы, Е — напряжённость электрического поля (В/м), Б — индукция электрического поля (Кл/м2), Ы — напряжённость магнитного поля (А/м), Б — индукция магнитного поля (Тл), ^ — плотность стороннего электрического тока (А/м2), р — плотность электрических зарядов (Кл/м3), ш — циклическая частота (Гц), е — диэлектрическая проницаемость (Ф/м), ^ — магнитная проницаемость (Гн/м), а — электрическая проводимость (См/м), г — мнимая единица.
Исключив напряжённость магнитного поля Ы, перейдем к векторному уравнению Гельмгольца относительно напряжённости электрического поля Е и сформулируем краевую задачу
го%-1го1Е) + к2Е = —гшЗв. (6)
Е х п
= 0,
да
где к2 = гша — ш2е, Е — комплексная векторная функция, П — область моделирования, дП — граница области, ^ — вектор плотности тока, ^у ^ = 0, так как в качестве источника возбуждения электромагнитного поля рассматривается соленоидальный контур.
На границах между подобластями Г^ с различными электрофизическими характеристиками должны выполняться условия непрерывности [1]
[п х Е]
[п(а + гше)Е]
0, (7)
0.
Сохранение заряда в областях с разрывными электрофизическими характеристиками обеспечивается выполнением условия непрерывности (8). Следовательно, вычислительные схемы как для решения систем уравнений Максвелла, так и для решения уравнения Гельмгольца должны конструироваться таким образом, чтобы условия (7), (8) выполнялись с заданной точностью. Векторный метод конечных элементов позволяет выполнить эти требования, если векторные базисные функции имеют второй и более полный порядок [4].
Г-
1 г]
1 г]
2. Векторная вариационная постановка
Пусть П — трёхмерная область, состоящая из неоднородных по физическим свойствам подобластей с липшицнепрерывной границей дП. Введём гильбертовы пространства векторных комплексных функций [5]
H(rot, П) = {v е [¿2(П)]3 : rot v е [¿2(П)]3},
Я01(П) = {ф е Ь2(П) : gradф е [£2(П)]3,ф
= 0
дП
H0(rot, П) = {v е H(rot, П) : v х n =0}
дП
со скалярным произведением и нормой [1].
Для задачи (6) сформулируем вариационную постановку: для Js е [Ь2(П)]3 найти E е H0(rot; П) такое, что для Vv е H0(rot; П) выполняется
(^-1rot E, rotv) + (k2E, v) = —¿(J v). (9)
Закон сохранения заряда имеет вид
div (Js + aE) = — гшр.
Так как div Js = 0 и р = div eE, то уравнение сохранения заряда будет следующим:
div (a + ¿^e)E = 0. (10)
В соответствии с диаграммой De Rham'a [6] имеет место свойство вложения grad ф е H0(rot; П), Vф е H0(grad; П). Возьмем v = grad ф е H0(rot; П). Тогда (9) имеет вид
(^-1 rot E, rot grad ф) + (k2E, grad ф) = — ¿(wJs, grad ф).
Поскольку rot grad ф = 0, div Js = 0, то в результате получим
((w2e + ¿wa)E, grad ф) = 0.
Следовательно, при выполнении условий вложения решение вариационной задачи удовлетворяет закону сохранения заряда (10) в слабом смысле.
3. Дискретный аналог вариационной задачи
Пусть задано разбиение области П на непересекающиеся тетраэдры, на каждом из которых определена иерархическая векторная базисная функция второго полного порядка. Тогда дискретный вариационный аналог в подпространстве Hh(rot; П) С H0(rot; П) задачи (9) формулируется в виде
Ou-1rot Eh, rot vh )п + (k2Eh, vh)n = — ¿(J vh)n.
Для дискретных подпространств Hh(rot; П) и Hh(grad; П) свойство вложения также будет выполняться:
vh е H0h(grad ; П) ^ grad uh е H0h(rot; П). Следовательно, закон сохранения заряда в этом случае выполняется в слабой форме.
Векторная функция E Е H^ (rot; П) может быть представлена как
Eh = ^ uiWi(x), (11)
ies
где щ — веса в разложении по базису Б — множество индексов степеней свободы, ж(ж1,ж2,ж3) — координаты в пространстве Я3.
Пространство Н^(гоЬ, П) является прямой суммой подпространств [7]:
Hh(rot; П) = N0h(rot; П) 0 (N0h(rot; П))х,
где Nq (rot; П) — ядро rot-оператора, (Nq (rot; П))х — его ортогональное дополнение. Для выполнения условий непрерывности (7) и (8) необходимо использовать базис второго типа, т. е. состоящий из роторных базисных функций, принадлежащих пространству (N^(rot; П))^ и обеспечивающих непрерывность тангенциальных компонент поля E (7), и градиентных базисных функций из пространства N^(rot;n), отвечающих за скачок нормальной компоненты поля (8) и выполнение закона сохранения заряда (10).
Таким образом, базисные функции в этом случае являются прямой суммой базисных функций двух типов: первый тип — из пространства (N^(rot; П))х, второй — из пространства N^(rot; П) [4, 8]. Иерархический базис второго полного порядка конструируется последовательным добавлением базисных функций более высокого порядка к роторным и градиентным базисным функциям первого порядка [9, 10]. Из дискретной вариационной формулировки получим СЛАУ
Au = f, (12)
где
Aij = (^-1rot wi, rot wj)n + (k2wi, wj)n, fi = -i(wJs, wi).
Структура матрицы A определяется, во-первых, векторными функциями, формирующими базисы полного порядка, во-вторых, тем, что базисы являются иерархическими по построению. Так как иерархический векторный базис состоит из функций первого и второго порядка, то структура матрицы A будет иметь вид
A
C11 C12 C21 С22
'13)
где Сц — блок матрицы, соответствующий дискретизации вариационной задачи (9) на базисе г-го порядка, С^ — блоки матрицы, соответствующие связи базисных функций г-го порядка с базисными функциями ]-го порядка при г = ]. Блоки С^ имеют блочную структуру, поскольку состоят из роторных и градиентных базисных функций
Cij
Du D12 D21 D22
'14)
где Dkk — блок, соответствующий дискретизации вариационной задачи (9) на базисе k-го типа, Dkl — блоки, соответствующие связи базисных функций k-го типа с базисными функциями 1-го типа (k = 1, 2; I = 1, 2).
Блочная структура матриц Cij (i = j) определяет связь базисов разного порядка и разных типов.
4. Методы декомпозиции с налеганием
Разобьём расчётную область О на т пересекающихся подобластей Ог
т
О = у О, V О = к :Ок П О = 0.
г=1
Область налегания подобласти к на подобласть 1 обозначим как Окг • Тогда Г кг — граница налегания подобласти к на подобласть 1, а Ггк — граница налегания подобласти 1 на подобласть к (рис. 1). В этом случае задача (6) эквивалентна двум подзадачам с условиями "сшивки" на интерфейсах Гкг и Ггк [11]
[п х Е] = 0, [п х Е]
= 0.
Данный алгоритм "сшивки" решения между подобластями соответствует схеме Дирихле — Дирихле для условия непрерывности тангенциальных компонент напряжённости электрического поля Е.
Мера налегания подобластей, нумерация геометрических примитивов (узлы, ребра, грани, тетраэдры) в каждой подобласти полностью определяют структуру СЛАУ (12) и число обусловленности матрицы.
г
1к
4.1. Особенности декомпозиции области
При построении геометрической декомпозиции области на подобласти необходима специальная процедура декомпозиции области моделирования, поэтому для каждого класса задач следует разрабатывать специальную процедуру разбиения области.
Рассмотрим три класса задач: 1) замкнутый токовый контур в однородной слабо проводящей области, 2) наземная электроразведка, 3) индукционный каротаж (сква-жинная электроразведка).
1. Для верификации алгоритма декомпозиции была решена задача в однородной области. Данная задача является модельной, так как не содержит особенностей, связанных с частичной непрерывностью решения на границе раздела сред с различными физическими свойствами.
2. Наземная электроразведка в слоистой среде (рис. 2), источник электромагнитного поля — петля, диапозон частот от 100 Гц до 100 кГц.
а--
Рис. 1. Пример разбиения области на подобласти
Рис. 2. Расчётная область: Г — генераторная катушка, 1, 2,3, 4 — горизонтальные слои, различные по физическим свойствам
3 2 1 а < к -Г
4 -И!
ЧИ2
Рис. 3. Скважина: 1 — зонд, 2 — скважина с буровым раствором, 3 — вмещающая среда, 4 — плохо проводящий слой, Г — генераторная катушка, И1, И — измерительные катушки
3. Индукционный скважинный каротаж (рис. 3), частота тока в источнике возбуждения поля 14 МГц. Многовариантность этой задачи связана с движением зонда внутри скважины.
4.2. Геометрическая декомпозиция области
При распределённом источнике поля возможно такое разбиение на подобласти, при котором матрицы в них будут иметь примерно одинаковую размерность. При таком разбиении СЛАУ (12) в каждой из подобластей решается примерно одинаковое время, что является наилучшим вариантом для параллелизации итерационного процесса.
В настоящей работе рассматривается локальный источник поля, что приводит к трём ограничениям:
1) источник поля должен находиться в одной подобласти и не касаться её границ. Это условие следует из схемы Дирихле — Дирихле, так как иначе при склейке решений правая часть исходного СЛАУ не будет соответствовать источнику возбуждения
в решаемой задаче;
2) граница области с источником должна отступать от него не менее, чем на 3-4 диаметра петли, поскольку поле в этой подобласти имеет наибольший градиент;
3) СЛАУ в подобластях, не содержащих источник, будут иметь нулевую правую часть (кроме граничных условий склейки), поэтому уточнение решения в них будет происходить быстрее, чем в области с источником.
Для автоматического разбиения области на подобласти с максимально возможным учётом накладываемых ограничений были разработаны следующие алгоритмы: а) декомпозиция на вложенных "сферах", б) декомпозиция по слоям. При декомпозиции с помощью сфер (рис. 4) сначала выбирается центральный тетраэдр (содержащий точку центра источника), затем постепенно наращиваются слои относительно предыдущих тетраэдров (алгоритм 1). В результате получим множество слоёв вокруг центрального тетраэдра. Из этих слоёв составляются подобласти с налеганием при условии максимальной сбалансированности подобластей по числу тетраэдров, что необходимо для обеспечения одинакового времени решения СЛАУ во всех подобластях. При декомпозиции по слоям в качестве начального слоя выбираются тетраэдры, узел, ребро или грань которых принадлежат одной из границ расчётной области. Дальнейшее построение слоёв тетраэдров такое же, как и в случае декомпозиции на "сферах".
Рис. 4. Декомпозиция на вложенных сферах
Алгоритм 1. Построение согласованных слоёв на вложенных сферах.
Th — список тетраэдров, n — текущее число тетраэдров. k — номер центрального тетраэдра.
Current_layer — множество тетраэдров в текущем слое. Next_layer — множество тетраэдров в следующем слое. Копировать Th[k] в Current_layer. Удалить Th[k]. n = n — 1.
Mas_layer — массив списков слоёв.
Повторять пока Th — не пустое. {
Для i = от 1 до n повторять {
Если Th[i] касается Current_layer, то {
Копировать Th[i] в Next_layer Удалить Th[i]
}
}
Копировать Current_layer в Mas_layer Current_layer = Next_layer
}
Next_layer = NULL.
Из полученных наборов слоёв тетраэдров составляются подобласти декомпозиции с учётом меры налегания (количество слоёв налегания), размерностей слоёв (балансировка матриц по размерности) и локальности источника.
Структура матрицы не зависит от способа геометрического разбиения области (сферы, слои), так как при обоих подходах подобласти пересекаются только с предыдущей и последующей подобластями, если такие имеются.
5. Структура матрицы
Рассмотрим алгебраическую декомпозицию, ассоциированную с матрицей А размерности п х п дискретного аналога вариационной задачи (9).
Структура матрицы будет иерархической, т. е. она содержит блоки, соответствующие подобластям декомпозиции (П,) [12], каждый из которых имеет вложенную блочную структуру, связанную с размерностью и типом базисных функций. Введём обозначения
^ i - П i
1 \ / т
и пА и и
\,+1
1,..., т.
п i р| П
1, ... , т,
3
1 , ... , т.
Рассмотрим блочную структуру матрицы А на примере трёх подобластей разбиения (к — 1,к,к +1) (рис. 5). Примем следующую последовательность нумерации базисных функций: сначала пронумеруем базисные функции, соответствующие подобласти Пк—1, затем — подобластям Пк—1к, Пк, Пкк+1, в последнюю очередь — подобласти П ^+1. Тогда матрица имеет вид
А
' Ак—1
А к— 1 Ак—1к
0 0 0
Ак— Ак— 1к 1 0 0 0
Ак— 1к лк Ак—1к 0 0
Ак— 1к Ак Акк+1 0
0 лк Акк+1 Акк+1 Ак+1 Акк+1
0 0 лкк+1 Ак+1 Ак+1 _
15)
где А, — блок, соответствующий дискретизации уравнения на подобласти П,, А^ — блок, соответствующий связи подобластей П , и П ^ и вычисляемый при дискретизации исходного уравнения в области П^. Свойства блока А^ полностью определяются мерой
ЭР
У
' I ' / /
I йк-г
/
I
Т —'
\ ' \ '
/ / __
к+1
к + 1
Рис. 5. Разбиение расчётной области на три подобласти
Рис. 6. Блочная структура матрицы
налегания между подобластями. Лг1 и Лгк1 — блоки, соответствующие связи подобласти Пг с налеганием её на другие подобласти.
Каждой из подобластей в матрице Л соответствует свой макроблок Вг
Bi
Л Ai A12
Ai2
A
12
Bk
Ak-lk лк—1к Ak
0
Ak Ak—lk
Ak Ak Akk+1
0
Akk+1
A
kk+1
Bm
Am-lm A^-lm
Am
m— lm A
Из структуры блоков Bk следует, что блоки k — 1 и k + 1 могут обрабатываться независимо друг от друга. Каждый из блоков (15) в свою очередь в силу иерархичности базисов разного порядка (11) тоже будет иметь блочную структуру. В результате блочная структура матрицы A может быть представлена схемой, изображённой на рис. 6. Такая сложная блочная структура матрицы позволяет конструировать не только алгебраические методы декомпозиции области, но и алгебраические многосеточные методы [13] с выбором в качестве грубых подпространств пространства, натянутые на базисные функции низкого порядка и типа (уточнение решения в различных блоках по порядкам и типам).
Алгебраический метод декомпозиции определяется процедурой проектирования решения на подпространство меньшей размерности, в качестве которого может быть использовано пространство, соответствующее подобласти, или подпространство, натянутое на базис более низкого уровня. Для построения этого алгоритма введём проектор из всей области на подобласть
P : Я^(rot; Q) ^ H0h(rot; Q-).
Матричным аналогом проектора P будет прямоугольная матрица, состоящая из нулей и единиц, в каждом столбце и строке которой не более одной единицы. С помощью матрицы R размерности n х n определим n х n-матрицу Bj как сужение матрицы A на подобласти Q
B,- = r-ART . (16)
Матрицы Bj вычисляются дискретизацией уравнения (9) на i-й подобласти.
6. Мультипликативный метод Шварца
Мультипликативный метода Шварца для решения СЛАУ (12) представляет собой итерационный процесс с т подшагами для нахождения очередного приближения , если задано текущее приближение и-7 [12, 14],
В^ = Д 7 - , (17)
т = т + ¿г, г = 1,...,т, г — номер подобласти, (18)
где Вг имеет вид (16), гг — поправка решения на г-й подобласти. Один полный шаг метода Шварца состоит из т подшагов (по числу подобластей), поэтому решение т имеет дробный индекс: целая часть ] — номер текущего полного шага, дробная часть г
--уточнённое решение на текущем полном шаге в г-й подобласти.
т
Определим матрицу перехода для одного полного шага. Так как
>-1
В-1 Д 7 -
то из (18) следует
—т = —^ + ДгтВг-1Дг 7 - А—
т = (I - ДТВг-1ДгА)и7+^ + Вг-1Дг/.
и
Введём обозначения: Сг = (I - Вг 1ДгА) — матрица перехода для подшагов, дг Вг-1Дг/ — вклад от правой части. В этом случае получим
т = С—^ + дг. (19)
Для следующего подшага (подобласти) (19) будет иметь вид
7+ ^ + 1 .-у 7 +
т = Сг+1М-+т + дг+1,
^ = Сг+1Сги-+^ + Сг+1дг + дг+1 = Сг+^г-^ + $+1. Выполнив т подшагов, получим
и-+1 = СтСт-1...С1М7 + д-+1.
Иными словами, матрица перехода между двумя полными шагами является произведением матриц перехода по подшагам [15], поэтому данный метод называется мультипликативным. При алгебраическом представлении методов декомпозиции области технологическая (программная) реализация одинакова как для функциональной, так и для геометрической декомпозиции. Основная сложность — построение операторов проектирования и сужения матрицы исходной СЛАУ. В случае функциональной декомпозиции проектор полностью определяется базисами подпространства и исходного пространства. При построении геометрической декомпозиции с налеганием исходную расчётную область необходимо разбить на подобласти.
Параллельная реализация мультипликативного метода Шварца. Методы декомпозиции области имеют естественную параллельную структуру (подобласти, не налегающие друг на друга, обрабатываются параллельно). Построенные методы автоматической декомпозиции при наличии чётного числа подобластей позволяют строить параллельную версию мультипликативного метода.
При декомпозиции области на подобласти предложенными методами имеют место следующие сценарии:
1) первая подобласть пересекается только со второй подобластью;
2) последняя подобласть пересекается только с предыдущей подобластью;
3) внутренние подобласти пересекаются со следующей и предыдущей подобластями.
Разобьём все множество 2п подобластей на п пар (рис. 7): г-я пара состоит из 2г — 1
и 2 г подобластей. Каждая из п пар соответствует одному узлу кластера. Подобласти с нечётными номерами могут обрабатываться независимо друг от друга, так как не имеют пересечений. Аналогично с чётными номерами. Внутри одной пары подобласти обрабатываются последовательно.
В этом случае алгоритм параллельной реализации на г-м узле решения уравнения Аи = f с точностью 7 имеет следующий вид.
Алгоритм 2 (декомпозиции области на одном узле кластера).
1. Выбираем начальное приближение и0г—1 = 0, и^ = 0, г = 1,..., п — номер пары
подобластей.
2. г
2г—1 0
= у 2г—1 — ^2
г—1„2г—1 г2г ип , I п
= у2г _ А2г
и
2г
3. Пока ,, „,, > 7
повторять.
4.1. г2г—1 = £(А2г
4.2. и2г—1 = и2г—1 + г
—1 г2г—1
2г 1
72 ).
4.3. < = < + Д2гЯТг—1 г
4.4. Если г > 1, то г2г—1
2г 1
передаем на г—1-й.
4.1. Если г < п, то принимаем данные от г + 1-го узла: г2г+1.
4.2. и0г = < + ^£+1
2г+1
5.1. г2г = £(А2г,г2г,72).
2г
5.2. и0г = и0г + г
5.3. и00г—1 = и0г—1 + Я*—1ЯТ г2г.
5.4. Если г < п, то г2г передаем на г + 1-й.
5.1. Если г > 1, то принимаем данные
2г 2
от г + 1-го узла: г
5 2 и2г—1 = и2г—1 + О оТ „2г—2 5.2. и0 = и0 + Л2г—1 Л2г_2г .
Здесь О — оператор проектирования на ]-ю подобласть.
0
и
Рис. 7. Пример декомпозиции для параллельной реализации; 1,2 — подобласти
Таким образом, алгоритм уточнения решения на одном узле кластера состоит из двух последовательных шагов 4 и 5. Каждый из них состоит из двух параллельных частей: уточнение решения в одной подобласти и обмен данными с соседним узлом. S(А, г, 7) — итерационный решатель с матрицей СЛАУ А, вектором правой части г
и точностью 7 — ^к_1 < £ << 1. Уточнение решения в подобластях выполняется
с небольшой точностью. Так как гс^-оператор имеет ядро большой размерности, то матрица дискретизации будет знаконеопределённой, следовательно, стандартные пред-обусловленные итерационные методы либо не будут сходиться, либо будут сходиться чрезвычайно медленно. Поэтому в качестве решения СЛАУ на подобласти используется двухуровневый итерационный решатель с коррекцией на пространстве ядра гс^опера-тора [16].
7. Алгебраический многосеточный метод
При локальном источнике поля существует проблема быстрой передачи информации между подобластями, так как ненулевая правая часть (источник электромагнитного поля) сосредоточена в одной из подобластей. Для её решения разработан алгебраический многосеточный алгоритм с использованием в качестве грубых пространств базиса более низкого уровня иерархии. Построение такого многосеточного алгоритма обеспечивается структурой матрицы, блоки которой определяются порядками и типами базисных функций (13), (14). Введём проекторы
Ргд : №22 ^ Wi3,
где — пространство, натянутое на базис ¿-го порядка 3-го типа. Матрицы, соответствующие данным проекторам, обозначим как Я^. Эти матрицы состоят из нулей и единиц и имеют блочную структуру
Яг] — (Е 0) ,
где Е — единичная матрица размерности &ш(Е) — ), 0 — нулевая матрица.
Размерность проекционной матрицы ) х ^ш(Ш22)).
Алгоритм 3. Алгебраический многосеточный метод.
1. и0 — начальное приближение, 7 — точность решения СЛАУ.
2. го — f — Ащ.
3. Для г — 1, 2,....
4. дп — Япг,_1.
5. ¿п — 5(Ап,дп,71).
6. и,_3 — и,_1 + Я^ц.
7. г,_ з — f — Аи,_ з.
, 4 * , 4
8. д12 — Я12Г, _ з.
9. ¿12 — 5(А12,д12,72).
10. щ_2 — щ_з + Я^2г12.
11. г,_ 2 — f — Аи,_ 4.
12. д21 — Я21г,_ 2.
13. ¿21 — 5(А21, д21,73).
14. и,_ 1 — и,_ 2 + Я^^ъ
15. г, 1 — / — Аи, 1.
' 4 1 4
16. д22 — г,_4.
17. ¿22 — 5(А, д22,74).
18. и, — и,_ 1 + ¿22.
' 4
19. г, — f — Аи,.
20. Если ||г,|| < 7||f ||, то вернуть и,, иначе — увеличить г.
Здесь Б(А^, д, 7) — решение СЛАУ на базисе ¿-го порядка 3-го типа с матрицей А^ — Я^ АЯТ, вектором правой части д и заданной величиной 7. В настоящей работе в качестве Б (...) используется мультипликативный метод Шварца.
Таким образом, внешним решателем является алгебраический многосеточный метод, на каждой подытерации которого используется мультипликативный метод Шварца. Данный метод в свою очередь на каждой из подобластей уточняет решение двухуровневым итерационным методом с коррекцией на пространстве ядра гс^оператора. В результате общая структура построенного алгоритма решения СЛАУ имеет следующий вид.
1. Многосеточный метод на базисах разного порядка и разного типа.
1.1. Решение СЛАУ на базисе первого порядка первого типа методом декомпозиции области.
1.1л. Решение СЛАУ на г подобласти двухуровневым итерационным решателем.
1.2. Решение СЛАУ на базисе первого порядка второго типа.
1.2л. Решение СЛАУ на г подобласти двухуровневым итерационным решателем.
1.3. Решение СЛАУ на базисе второго порядка первого типа.
1.3л. Решение СЛАУ на г подобласти двухуровневым итерационным решателем.
1.4. Решение СЛАУ на базисе второго порядка второго типа.
1.4л. Решение СЛАУ на г подобласти двухуровневым итерационным решателем.
8. Вычислительный эксперимент
а. Моделирование электрического поля в слоистой среде с локальным наземным источником поля (петля). Разрез расчётной области представлен на рис. 2.
Радиус генераторной катушки г — 1 м. Геометрические параметры слоёв (высота слоя): к1 — 20, к2 — 100, к3 — 40, к4 — 20 м. Сила тока I — 1 А, частота источника поля f — 1 кГц. Физические параметры среды: £, — £0, — ^0, I — 1, 2, 3, 4, а1 — 0.1, а2 — 1, а3 — 0.2, а4 — 0.001 Сим/м.
Декомпозиция на сферы.
Размерности СЛАУ:
Без декомпозиции: 4059945.
Декомпозиция на две подобласти с налеганием — один слой тетраэдров: 2302346, 2016234.
Декомпозиция на две подобласти с налеганием — два слоя тетраэдров: 2819082, 2016234.
Декомпозиция на две подобласти с налеганием — три слоя тетраэдров: 2819082, 2359383.
Декомпозиция на четыре подобласти с налеганием — один слой тетраэдров: 1780701, 1614330, 1290852, 1063206.
Декомпозиция на четыре подобласти с налеганием — два слоя тетраэдров: 1780071, 2203722, 1880247, 1063206.
Декомпозиция на четыре подобласти с налеганием — три слоя тетраэдров: 2135829, 1837413, 1416888, 1063206.
Декомпозиция на восемь подобластей с налеганием — один слой тетраэдров: 589968, 615849, 646362, 931077, 1118520, 772815, 584070, 631113.
Декомпозиция на восемь подобластей с налеганием — два слоя тетраэдров: 725778, 615849, 795801, 1125330, 1397976, 1156746, 827554, 758229.
Декомпозиция на восемь подобластей с налеганием — три слоя тетраэдров: 288585, 485673, 790731, 1223742, 1808583, 1634001, 939411, 758229.
Время решения СЛАУ (последовательная реализация), с
Число Налегание между подобластей подобластями 1 2 3
2 2053 1681 2906
4 3028 2013 2504
8 7255 5199 4613
Время решения СЛАУ (параллельная реализация), с
Число Налегание между подобластей подобластями 1 2 3
4 2350 1934 1878
8 2234 2135 2973
Декомпозиция на слои. Размерности СЛАУ:
Декомпозиция на две подобласти с налеганием — один слой тетраэдров: 2107113, 2134371.
Декомпозиция на две подобласти с налеганием — два слоя тетраэдров: 2292501, 2134371.
Декомпозиция на две подобласти с налеганием — три слоя тетраэдров: 2299501, 2266518.
Декомпозиция на четыре подобласти с налеганием — один слой тетраэдров: 1074186, 1198812, 1130286, 1323507.
Декомпозиция на четыре подобласти с налеганием — два слоя тетраэдров: 1132574, 1378584, 988722, 1592232.
Декомпозиция на четыре подобласти с налеганием — три слоя тетраэдров: 1485225, 1510842, 1357296, 1323507.
Декомпозиция на восемь подобластей с налеганием — один слой тетраэдров: 697833, 791907, 829056, 873480, 1007166, 1040385, 542997, 1252122.
Декомпозиция на восемь подобластей с налеганием — два слоя тетраэдров: 725778, 615849, 795801, 1125330, 1397976, 1156746, 827554, 758229.
Время решения СЛАУ (последовательная реализация), с
Число Налегание между подобластей подобластями 1 2 3
2 3082 2462 1758
4 5343 3978 3502
8 9715 6319
Время решения СЛАУ (параллельная реализация), с
Число Налегание между подобластей подобластями 1 2 3
4 3076 2751 2603
8 3354 2535
"Ручная" декомпозиция.
Размерности СЛАУ:
Без декомпозиции: 5544136.
Декомпозиция на четыре подобласти: 1643649, 1292025, 2810763, 1018665.
Время решения СЛАУ, с
Число Реализация
подобластей последовательная параллельная
1 1277 —
4 1574 952
При построении геометрической декомпозиции необходимо выполнение следующего условия: мера налегания должна быть согласована с длиной волны, что необходимо для сохранения физического смысла распределения вычисленного векторного поля. Поэтому налегание в один слой тетраэдров не всегда даёт хороший результат. При использовании декомпозиции на слои при такой внутренней геометрии расчётной области получается лучший результат, чем в случае декомпозиции на вложенных "сферах", так как больше соответствует структуре области моделирования. Однако построение автоматической декомпозиции не даёт выигрыша во времени по сравнению с решением этой задачи с помощью "ручной" декомпозиции: границы области являются липшицнепре-рывными, матрицы максимально сбалансированы по подобластям. Поэтому в данном случае получается выигрыш во времени. Но такой подход является трудоёмким и при усложнении внутренней геометрии области моделирования подобрать размерности подобластей и меры налеганий становится всё труднее.
б. Скважина, в которой находится зонд. Частота источника тока ш =14 МГц, сила тока I =1 А. Расчётная область представлена на рис. 3. Численное моделирование проводилось при = е0, = для всех вариантов.
Номер области 1 2 3 4 аг 0 5 0.1 0.01
Декомпозиция на сферы.
Размерности СЛАУ:
Без декомпозиции: 2956671.
Декомпозиция на две подобласти с налеганием — один слой тетраэдров 1274712.
Декомпозиция на две подобласти с налеганием — два слоя тетраэдров 1616706.
Декомпозиция на две подобласти с налеганием — три слоя тетраэдров 1616706.
: 2106912, : 2106912, : 2400201,
Декомпозиция на четыре подобласти с налеганием — один слой тетраэдров: 1428849, 1063317, 943101, 628572.
Декомпозиция на четыре подобласти с налеганием — два слоя тетраэдров: 1428849, 1314372, 1285095, 934752.
Декомпозиция на четыре подобласти с налеганием — три слоя тетраэдров: 1116801, 1696656, 1716660, 934752.
Декомпозиция на восемь подобластей с налеганием — один слой тетраэдров: 632067, 706545, 636309, 730290, 766947, 718242, 603141, 628572.
Декомпозиция на восемь подобластей с налеганием — два слоя тетраэдров: 236826, 332079, 636759, 1018593, 1314372, 1285095, 810450, 381861.
Декомпозиция на восемь подобластей с налеганием — три слоя тетраэдров: 64398, 160590, 391128, 955374, 1696656,1716660, 873222, 223779.
Время решения СЛАУ (последовательная реализация), с
Число Налегание между подобластей подобластями 1 2 3
2 2851 1978 1569
4 3622 2964 2941
8 6440 4303 -
Время решения СЛАУ (параллельная реализация), с
Число Налегание между подобластей подобластями 1 2 3
4 3054 1762 2387
8 2778 2451 -
"Ручная" декомпозиция представлена на рис. 8.
Размерности СЛАУ:
Без декомпозиции: 3017698.
Декомпозиция на две подобласти: 1497276, 1876451.
Время решения СЛАУ, с
Число Время
подобластей решения
1 862 2 723
В данном случае нет ярко выраженной слоистой структуры, поэтому для автоматической декомпозиции используется только декомпозиция на вложенные сферы. Декомпозицию на восемь подобластей с налеганием в три слоя тераэдров здесь построить нельзя, поскольку общее число слоёв тетраэдров не достаточно. Как и в слоистой среде, увеличение меры налегания улучшает процесс сходимости, но при этом ухудшает параллельные свойства алгоритма (увеличивает объём обмениваемыми данными). Построение "ручной" декомпозиции для такой расчётной области сложнее, чем в слоистой
Рис. 8. Декомпозиция на две подобласти
среде, поэтому декомпозиция строится только на две подобласти. Но и в этом случае при правильно построенных налеганиях и сбалансированных матрицах подобластей получается выигрыш во времени, так как матрицы подобластей обрабатываются быстрее, чем одна большая матрица, соответствующая всей расчётной области.
Список литературы
[1] Баландин М.Ю., Шурина Э.П. Векторный метод конечных элементов: Уч. пособие. Новосибирск: Изд-во НГТУ, 2001.
[2] Rodrigue G., White D. A vector finite element time-domain method for solving Maxwell's equations on unstructured hexahedral grids // SIAM J. Sci. Comput. 2001. Vol. 35. P. 315-341.
[3] Нечаев О.В., Шурина Э.П. Многосеточный алгоритм решения векторным методом конечных элементов трёхмерного уравнения Гельмгольца // Матем. моделирование. 2005. Т. 17, № 6. C. 92-102.
[4] Webb J.P. Edge elements and what they cando for you // IEEE Trans. on Magnetic. 1993. No. 2. P. 1460-1465.
[5] Nedelec J.C. Mixed finite elements in R3 // Numer. Math. 1980. No. 3. P. 315-341.
[6] Schwarzbach C. Stability of Finite Element Solutions to Maxwell's Equations in Frequency Domain: Dis. Dr. Rer. Nat. Gorlitz, 2009.
[7] Hiptmair R. Multigrid methods for Maxwell's equations // SIAM J. Nymer. Anal. 1998. No. 1. P. 204-225.
[8] Bossavit A. Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements. San Diego, USA: Acad. Press Ltd., 1998.
[9] Lars S. Andersen, John L. Hierarchical tangential vector finite elements for tetrahedra // IEEE Microwave and Guide Wave Lett. 1998. No 3. P. 8.
[10] Nedelec J.C. A new family of mixed finite elements in R3 // Numer. Math. 1986. No. 50. P. 97-81.
[11] Васильевский Ю.В., ОльшАнский М.А. Краткий курс по многосеточным методам и методам декомпозиции области. М.: МАКС Пресс, 2007. 104 c.
[12] Gander M.J. Schwarz methods over the course of time // Electronic Trans. on Numer. Anal. 2008. Vol. 31. P. 228-255.
[13] Mense C., Nabben R. On algebraic multilevel methods for non-symmetric systems-convergence results // Ibid. 2008. Vol. 300. P. 323-345.
[14] Toselli A., Widlung O. Domain Decomposition Methods-Algorims and Theory. Springer, 2004.
[15] Benzi M., Frommer A., Nabben R., B.Szyld D. Algebraic theory of multiplicative Schwarz methods // Numer. Math. 2001. No. 89. P. 605-639.
[16] Nechaev 0., Shurina E., Botchev M. Multilevel iterative solvers for the edge finite element solution of the 3D Maxwell equation // Comput. and Math. With Appl. 2008. No. 10. P. 16.
nocmynujia e peda'^uw 30 non6pn 2012 s.
