CC BY
8Научный вестник НГТУ. - 2011. - № 3(44)
УДК 519.63
Конечноэлементное моделирование трехмерных магнитотеллурических полей с применением технологии деревьев-кодеревьев и постановки с совместным использованием векторного и скалярного потенциалов*
П.А. ДОМНИКОВ, С.В. КИРЕЕВА, М.Г. ПЕРСОВА, Ю.Г. СОЛОВЕЙЧИК
В статье рассматривается моделирование трехмерных магнитотеллурических полей с использованием векторного метода конечных элементов. Проводится сравнение эффективности вычислительной схемы, основанной на использовании конечноэлементной постановки для вектор-потенциала A с исключением нуль-ядра rot rot-оператора методом tree-cotree, и схемы с применением A-V-постановки, в которой для аппроксимации вектор-потенциала A используются реберные базисные вектор-функции, а для скалярного потенциала V — скалярные узловые базисные функции.
Ключевые слова: векторный метод конечных элементов, метод деревьев-кодеревьев, магнитотеллурическое зондирование.
ВВЕДЕНИЕ
Выполнение работ по проектированию и интерпретации данных геоэлектромагнитных исследований в реальных трехмерных средах требует проведения вычислительно трудоемких расчетов, причем с повышением требований к качеству выполнения электроразведочных работ растут требования к эффективности расчетов трехмерных геоэлектромагнитных полей.
Одним из наиболее эффективных подходов к уменьшению вычислительных затрат при решении трехмерных задач геоэлектромагнетизма является выделение основной части поля горизонтально-слоистой среды [1]. Однако вычислительные затраты, требуемые для расчета добавочного трехмерного поля в сложных средах, содержащих большое количество трехмерных включений, могут оставаться существенными.
Среди задач геоэлектроразведки заметное место занимают магнитотеллурические зондирования (МТЗ) [2-4], требующие расчетов гармонических по времени электромагнитных полей в частотном диапазоне примерно от 10_4 до 500 Гц. Ранее были проведены исследования [5, 6], которые показали, что при использовании узлового и векторного метода конечных элементов (МКЭ) время расчетов трехмерных магнитотеллурических полей существенно зависит от частоты источника. Причем основные вычислительные затраты приходятся на несколько задач, соответствующих частотам из среднего диапазона, как правило, порядка 0.1 - 1 Гц. Это обстоятельство также отрицательно влияет на эффективность распараллеливания задачи МТЗ по гармоникам [6].
В данной статье приводятся результаты исследований, направленных на сокращение вычислительных затрат при моделировании трехмерных магнитотеллурических полей с использованием векторного МКЭ.
* Получена 26 апреля 2011 г.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг.
1. МАТЕМАТИЧЕСКИЕ МОДЕЛИ И ВАРИАЦИОННЫЕ ПОСТАНОВКИ
Математические модели для описания трехмерных магнитотеллурических полей, рассматриваемые в данной работе, основаны на использовании технологии выделения основной
части поля горизонтально-слоистой среды [1]. Вектор-потенциал магнитного поля A (являющийся комплекснозначной величиной) представляется в виде суммы A = Л" + Ла, где величина Л" (нормальное поле) описывает поле в горизонтально-слоистой среде, а величина Ла представляет собой добавочное поле от трехмерных неоднородностей. При этом вектор-
потенциал нормального поля Л" может быть найден из решения одномерной задачи
1 ^ Л (+ то" (2) Л" (2) = 0, (1)
М0 dz2
1 dAn (z)
цо dz
= 1 + 0 • i, (2)
z=0
где под An (z) понимается x - или y -компонента вектор-потенциала An (поскольку при решении трехмерной задачи МТЗ обычно используется два ортогональных направления сторонних токов, направленных вдоль осей x и y соответственно), Мо - магнитная проницаемость вакуума, cn (z) - удельная электрическая проводимость горизонтально-слоистой среды, ю = 2tcv - круговая частота, v - частота поля, i - мнимая единица.
Добавочное поле, описываемое вектор-потенциалом Aa, может быть найдено из решения следующего векторного уравнения
— rot rotAa + i®cAa = (о - on )En (3)
Цо
с однородными краевыми условиями первого рода (для касательных составляющих A), заданными на границе расчетной области. В (3) с = с( x, y, z) - функция, описывающая удельную проводимость среды с включенными в нее трехмерными неоднородностями, En = — iaAn - нормальная составляющая напряженности электрического поля во вмещающей горизонтально-слоистой среде. Индукция магнитного поля Ba и напряженность электрического поля E , соответствующие аномальному полю, определяются соотношениями:
Ba = rot Aa, Ea = —iaAa.
Эквивалентная вариационная постановка в форме Галеркина для уравнения (3) при условии, что пробная вектор-функция у принадлежит пространству Но (rot, Q) функций, квадратично интегрируемых вместе со своим ротором в области Q и имеющих нулевую касательную составляющую на ее границе SQ , имеет следующий вид:
— J rot Aa • rotydft + ia J oAa • ydft = J (о — on)En • ydft. (4)
Ц0 Q Q Q
В данной работе будут рассматриваться конечноэлементные аппроксимации на векторных параллелепипедальных элементах. Локальные базисные вектор-функции, определенные на шаблонном элементе [—1,1] х [—1,1]х [—1,1], имеют вид
\ = (Ф1 (У)Ф1 (z), 0, 0), у = (Ф2 (У)Ф1 (z), 0, 0), \уз = (ф1 (у)ф2 (z), 0, 0),
У4 = (Ф2 (У)Ф2 (z), 0, 0), \5 = (0, ф1 (х)ф1 (z), 0), у = (0, ф2 (х)ф1 (z), 0), (5)
/V /V /V (5)
\\7 = (0, Ф1 (Х)Ф2 (z), 0), = (0, Ф2 (Х)Ф2 (z), 0), = (0, 0, Ф1 (Х)Ф1 (y)), \\10 = (0, 0, Ф2 (х)Ф1 (y)), \1 = (0, 0, Ф1 (х)Ф2 (y)), \\12 = (0, 0, Ф2 (х)Ф2 (y)),
1 _ ^ 1 + ^
где Ф1 = и Ф2 = - одномерные скалярные функции. Глобальные базисные
вектор-функции строятся сшивкой локальных базисных функций, принадлежащих одному и тому же ребру конечноэлементной сетки.
Конечноэлементная аппроксимация на основе вариационной постановки (4) приводит к системе линейных алгебраических уравнений (СЛАУ)
^^ = ^ , (6) где комплекснозначная матрица P состоит из элементов pjk = bjk + iCjk, в которых элементы матрицы жесткости B = {bjk } и элементы матрицы массы C = [cjk } вычисляются по формулам
bjk = — J rot y j ■ rot ykdQ, (7)
^0 Q
c
jk
= ю|суj -ykdQ . (8)
Численные эксперименты, проведенные в [5, 6], показывают, что, во-первых, вычислительные затраты, требуемые для решения системы (6) итерационными методами, существенно возрастают при уменьшении частоты поля ю , во-вторых, использование предобусловливате-лей, основанных на неполной факторизации матрицы СЛАУ для решения системы (6), является неэффективным.
В научных публикациях довольно часто обсуждаются проблемы, связанные с медленной скоростью сходимости итерационных методов и с предобусловливанием при решении задач моделирования электромагнитных полей в технических устройствах с использованием векторного МКЭ (ВМКЭ) [7-12]. При этом отмечается, что одной из возможных причин медленной сходимости итерационных методов при решении СЛАУ, получаемых в результате ВМКЭ-аппроксимации задач электромагнетизма, является наличие большого нуль-ядра у оператора rotrot (состоящего из функций вида grad ф), ВМКЭ-аналогом которого является матрица
жесткости B = {bjk } .
Действительно, если расчетная область содержит непроводящие подобласти (т. е. подобласти с нулевым значением в них коэффициента удельной электрической проводимости с),
то решением уравнения (3) кроме функции Aa является любая функция вида A' = Aa + grad ф, где ф - некоторая скалярная величина, принимающая ненулевое значение в подобласти с с = 0 (и равная нулю в подобласти с с Ф 0 ). При этом магнитная индукция будет вычислена корректно: rot A' = rot( Aa + grad ф) = rot Aa = Ba.
Для устранения нуль-ядра rot rot -оператора существует технология деревьев-кодеревьев (tree-cotree technique) [7-12], суть которой заключается в следующем. В конечноэлементной сетке строится остовное дерево, вершинами которого являются узлы конечноэлементной сетки, а ребрами, соответственно, ребра конечноэлементной сетки. Пусть ne - число ребер в сетке, nn - число узлов в сетке, Nt - множество ребер сетки, принадлежащих остовному дереву, Nc - множество ребер, не вошедших в дерево (т. е. формирующих кодерево). При этом dim(Nt) = nn -1, dim(Ne) = ne - nn +1. После построения дерева (и, соответственно, коде-рева) базисные вектор-функции y j , j е Nt заменяются градиентами от скалярных базисных функций фj , j = 1,..., nn — 1, связанных со всеми узлами сетки за исключением узла, из кото-
рого началось построение дерева. Поскольку rotgrad ф= 0, то при использовании базисных вектор-функций вида grad ф элементы матрицы жесткости (7) дают нулевые вклады в глобальную матрицу конечноэлементной СЛАУ. В среде с с = 0 элементы матрицы массы (8) также равны нулю. Поэтому неизвестные, связанные с ребрами дерева, находящимися в среде с с = 0, можно исключить.
Однако в ряде работ (например, [9-11]) отмечается, что во многих случаях использование технологии деревьев-кодеревьев не приводит к улучшению сходимости итерационных методов решения конечноэлементной СЛАУ (несмотря на то, что матрица СЛАУ становится невырожденной). В работе [9] при решении задачи моделирования высокочастотного электромагнитного поля в волноводе использовался совершенно противоположный подход, заключающийся в том, чтобы не исключать функции вида grad ф, а явно добавить их в базис. В нашем исследовании мы применим аналогичный подход для решения трехмерной задачи магнитотеллурического зондирования с выделением основной части поля. Вариационная постановка в этом случае будет иметь вид
— J rot Aa • rotydQ + /ш J o(Aa + grad Va) • ydQ = J (о - о")Ё" • ydQ, (9)
^0 Q Q Q
/Ш J o(Aa + grad Va) • grad фdQ = J (о - о")Ё" • grad фdQ . (10)
Q Q
Уравнение (10) было получено умножением (4) на пробную функцию grad ф, где
ф е Hq , где Hq - гильбертово пространство скалярных функций, имеющих в Q суммируемые с квадратом первые производные и на SQ равные нулю. Напряженность аномального электрического поля и индукция аномального магнитного поля в этом случае вычисляются по
формулам Ba = rot Aa , Ea = -ira(Aa + gradVa). При этом вектор-потенциал Aa , как и ранее, раскладывается по базисным вектор-функциям (5), а скалярный потенциал Va представляется в виде линейной комбинации скалярных базисных функций, связанных с узлами конеч-ноэлементной сетки. Локальные скалярные базисные функции на шаблонном элементе [-1,1] х [-1,1] х [-1,1] имеют вид
ф1 = ф1 (х)ф (У )ф1 (z), ф2 = ф2 (х)ф (y )ф1 (z),
фз = ф1 (х)ф2 (У)ф1 (z), ф4 = ф2(х)ф2(у)ф^), (Q1)
ф5 = ф1 (Х)ф1 (у)ф2 (Z), ф6 = ф2 (Х)ф1 (у)ф2 (Z), ф7 = ф1( Х)фг( у)фг( z), ф8 = ф2( Х)фг( у)фг( Z).
Следует учесть, что градиенты от скалярных базисных функций (11) являются линейными комбинациями базисных вектор-функций (5), например, grad <Pi = - у5 - и т.д. Поэтому уравнения, получаемые при конечноэлементной дискретизации уравнения (10), являются линейно зависимыми по отношению к системе уравнений, получаемых при дискретизации уравнения (9).
2. ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ
Проведем вычислительный эксперимент на следующей геоэлектрической модели. Расчетная область состоит из вмещающей четырехслойной среды, содержащей первый слой толщиной 1 км с удельным сопротивлением 1400 Ом-м, второй слой толщиной 27 км с удельным сопротивлением 20000 Ом^м, третий слой толщиной 57 км с удельным сопротивлением
3000 Ом-м и последний слой с удельным сопротивлением 100 Ом-м. В данную среду помещены трехмерные неоднородности, размеры и удельное сопротивление которых приведены в табл. 1.
Построенная конечноэлементная сетка содержит 131040 прямоугольных параллелепипедов, 138966 узлов и 408817 ребер. Фрагмент сечения конечноэлементной сетки в плоскости YZ показан на рис. 1.
Рис. 1. Фрагмент сечения сетки с 9 объектами в плоскости х = 0
Таблица 1
Размеры и удельное сопротивление трехмерных неоднородностей
№ п/п Размеры по х , км Размеры по у, км Размеры по г,км Удельное сопротивление, Ом-м
1 -220. .400 -10 ..10 -16...-12 2
2 -70. .50 -100 ..100 -31,5...-28 2
3 -70. .50 -100 ..100 -63,5...-31,5 40
4 -220.. -193 -16 ..16 -2,5...-1 1400
5 -193.. -170 -16 ..16 -1...0 890
6 -193.. -170 -16 ..16 -3,4...-1 400
7 -193.. -170 -16 ..16 -8...-3,4 1000
8 -252.. -220 -16 ..16 -8...-5 1
9 -252.. -193 16.. .40 -8...-5 1
Для данной задачи была проведена серия расчетов аномального поля в частотном диапазоне от 0.00034 Гц до 500 Гц (всего 22 частоты, сторонний ток направлен вдоль оси X) при использовании следующих ВМКЭ-постановок.
1. Использовалась вариационная постановка (4) для вектор-потенциала Аа без использования технологии деревьев-кодеревьев и без введения скалярного потенциала. Далее такой подход будем называть А -постановкой.
2. Использовалась вариационная постановка (9)—(10) для вектор-потенциала Аа с введением скалярного потенциала Vа . Далее такой подход будем называть А-^постановкой.
3. Для А -постановки использовалась технология деревьев-кодеревьев, причем в дерево включались только ребра, находящиеся в подобласти с с = 0. Неизвестные, соответствующие ребрам дерева, исключались.
4. Для А -постановки использовалась технология деревьев-кодеревьев, причем в дерево были включены ребра, находящиеся как в подобласти с с = 0, так и в подобластях с с Ф 0 . Неизвестные, связанные с ребрами дерева, исключались, а в узлах сетки, лежащих в подобласти с с Ф 0 , был введен скалярный потенциал Vа .
Остовное дерево в расчетной области строилось в два этапа [8]. На первом этапе выбирался узел, лежащий на границе области, из которого начиналось построение дерева, и в дерево включались только ребра, принадлежащие границе области. На втором этапе после того как были пройдены все узлы, принадлежащие границе области, в дерево включались ребра, лежащие внутри расчетной области.
Аналогично при построении дерева только в подобласти с с = 0 сначала в дерево включались только ребра, лежащие на границе данной подобласти, а затем ребра, лежащие внутри нее. После этого ребра, лежащие на границе раздела сред с с = 0 и сФ0, исключались из дерева.
Так как остовное дерево не является единственным, были проведены вычислительные эксперименты для различных вариантов его построения. В первом случае для выбора вершин и ребер остовного дерева использовался обход в глубину. В этом случае в дерево оказывались включенными ребра, часть из которых направлена вдоль оси X, часть - вдоль оси У, часть -вдоль оси Z. Такое остовное дерево будем далее называть ХУ7-деревом. Поскольку на границе расчетной области заданы однородные краевые условия первого рода, для построения дерева по границе области во всех случаях использовался один и тот же алгоритм - обход в глубину. Отдельно были рассмотрены варианты построения остовного дерева, в которых внутри области в дерево включались ребра, направленные только вдоль одной из координатных осей: либо X, либо У, либо Z (такие остовные деревья далее будем называть Х-дерево, У-дерево и Z-дерево соответственно).
Число комплекснозначных неизвестных и заполненность матриц конечноэлементных СЛАУ при использовании различных постановок и вариантов построения дерева приведены в табл. 2.
Использование различных вариантов построения дерева приводит к конечноэлементным матрицам с разной заполненностью ненулевыми элементами, что влияет на время операции матрично-векторного умножения в итерационных алгоритмах решения СЛАУ. Поэтому при сравнении эффективности различных постановок будем приводить время решения СЛАУ.
На рис. 2-4 приведены графики зависимости времени решения конечноэлементных СЛАУ в зависимости от частоты поля V . В качестве метода решения СЛАУ использовался метод COCR [13] со сглаживанием невязки [14]. В качестве предобусловливателя использова-
т
лась либо главная диагональ (состоящая из комплексных чисел), либо неполное -разложение матрицы СЛАУ (элементы нижнетреугольной матрицы Ь также являются комплексными числами).
Критерием останова итерационного процесса служило уменьшение относительной нормы невязки в 105 раз или превышение 150 000 итераций (аварийный выход). При выходе по
Таблица 2
Размеры и заполненность матриц конечноэлементных СЛАУ
Число неизвестных Число ненулевых элементов в нижнем треугольнике матрицы
А-постановка 408 817 6 336 552
А^-постановка 502 467 12 350 211
XYZ-дерево (во всей области) 363 501 7 018 524
Х-дерево (во всей области) 363 501 7 013 680
Y-дерево (во всей области) 363 501 7 022 701
Z-дерево (во всей области) 363 501 7 056 083
XYZ-дерево (в подобласти) 363 501 5 136 067
Х-дерево (в подобласти) 363 501 5 134 609
Y-дерево (в подобласти) 363 501 5 139 037
Z-дерево (в подобласти) 363 501 5 154 096
д
А Й." Ф \
\ »
н./о. А
а, А ® \г д
V «Л V
л
с, о-е-о-о
•О— ХУг-1тее, СОСЩР^ - в - ХУг^гее, СОСК(ЬЬТ) ....д.... х-1тее СОСЩР^) —А— Х-1тее, СОСК(ЬЬТ) ....□•... у-1тее, СОСЩР^) —В— У^гее, СОСК(ЬЬТ) —О-.- г-1гее, СОСЩР^) ~0— 1-Пее, СОСК(ЬЬТ)
о о о о
1Л ГП ГМ <£> СО ГЧ СО Ш *Н *Н ~ тН гП Ш
V, Гц
Рис. 2. Время решения конечноэлементных СЛАУ в зависимости от частоты поля при различных вариантах построения дерева только в подобласти с ст = 0
критерию уменьшения относительной невязки выполнялся контроль точности получаемых в приемниках значений электромагнитного поля. Для всех сравниваемых постановок эта точность была примерно одинаковой.
Проанализируем результаты, полученные при использовании технологии деревьев-кодеревьев при построении остовного дерева только в подобласти с ст = 0 (рис. 3). Во-первых, время решения СЛАУ существенно зависит от выбора ребер, входящих в остовное дерево. При этом наилучший результат наблюдался при включении в остовное дерево ребер, направленных вдоль оси Z, а наихудший - при включении в остовное дерево ребер, направленных вдоль оси X (совпадающих с направлением стороннего тока). Во-вторых, использование неполной факторизации в качестве предобусловливания не приводило (за редким исключением) к уменьшению времени решения конечноэлементных СЛАУ, несмотря на исключение нуль-ядра. В третьих, с уменьшением частоты поля время решения СЛАУ очень быстро увеличивалось (с нескольких секунд на частоте 500 Гц до многих часов на самых низких частотах). При
этом решения задач для частот ниже 0.02 Гц удалось вычислить только с построением Г-дерева, 2-дерева и Л^-дерева и с использованием диагонального предобусловливания - в других случаях не удавалось уменьшить норму относительной невязки в 105 раз за 150 000 итераций (что соответствовало примерно 10 часам счета).
I, с
32768 16384 8192 4096 2048 1024 512 256 128 64 32 16 8 4 2 1
о ••"•■•о. □
'А. О.
Д. Ч. 'А О.
л-'' .-О'' О А'\
О .о И'' % *
»■•в
В N
о.
о. О-^-о-п
О. " V
О л , л
V-©--о
—О-- ХУг-иее, ....д.... х-иее, СОСК(01а§) ....□•... у-1тее, СОСК(Р1а§) —В— У-Шее, СОСЩШ") -•О-— Цхее, С0СЩ01а§) —0— Итее, СОСЩШ")
Ш СО СО
*н т-ч о
(Г) Ш Г*1
V, Гц
Рис. 3. Время решения конечноэлементных СЛАУ в зависимости от частоты поля при различных вариантах построения дерева во всей расчетной области
Рис. 4. Время решения конечноэлементных СЛАУ в зависимости от частоты поля при использовании А-постановки, А-К-постановки и построении Г-дерева во всей расчетной области
и 2-дерева в подобласти с ст = 0
Теперь проанализируем, как влияет на время решения задачи применение подхода, основанного на построении остовного дерева во всей расчетной области. При использовании данного подхода для частот ниже 0.6 Гц получить с требуемой точностью решения соответ-
ствующих конечноэлементных СЛАУ удалось только при построении Z-дерева (с использованием диагонального предобусловливания и LL -предобусловливания) и 7-дерева (с использованием диагонального предобусловливания). В других случаях не удавалось уменьшить
норму относительной невязки в 105 раз за 150 000 итераций. В частотном диапазоне от 0.01 Гц до 500 Гц наблюдался рост числа итераций при понижении частоты поля. Однако для частот в диапазоне 0.00034 - 0.01 Гц наблюдалось уменьшение времени решения СЛАУ с понижением частоты, и в данном частотном диапазоне время решения СЛАУ было значительно меньшим, чем при использовании подхода с построением остовного дерева только в подобласти с с = 0. Так, время решения СЛАУ для частоты 0.00068 Гц при построении 7-дерева только в подобласти с с = 0 составило 5 ч 11 мин (108670 итераций), а при построении 7-дерева во всей расчетной области - только 28 мин 37 с (7876 итераций).
Время, затраченное на решение СЛАУ при использовании A -постановки, A-V-постановки и наилучших результатов, полученных при использовании технологии деревьев-кодеревьев (построение дерева с включением ребер, направленных вдоль оси 7 во всей расчетной области и вдоль оси Z в подобласти с с = 0), приведено в виде графиков отдельно на рис. 4. В верхней части частотного диапазона 0.3-500 Гц использование A-постановки и подхода, основанного на построении остовного дерева в подобласти с с = 0, дает сравнимый результат. При этом подходы, основанные на использовании A-V-постановки и построении остовного дерева во всей расчетной области требуют большего времени решения СЛАУ. В нижней части частотного диапазона 0.00034 - 0.03 Гц наименьшее время решения конечноэлементных СЛАУ было получено при использовании A-V-постановки. При этом для частот ниже 0.005 Гц время решения СЛАУ для A-V-постановки с уменьшением частоты не увеличивалось. На самых низких частотах 0.00034 - 0.002 Гц применение подхода с построением остовного дерева во всей расчетной области позволило решить конечноэлементные СЛАУ за меньшее время, чем при использовании A-постановки, однако использование A-V-постановки дало еще лучший результат. Так, время решения СЛАУ методом COCR c диагональным предобусловливанием для частоты 0.00034 Гц при построении 7-дерева во всей расчетной области составило 19 мин 37 с (5448 итераций), при использовании A-постановки - 2 ч 36 мин (44540 итераций), при использовании A-V-постановки - 10 мин 54 с (1880 итераций).
ЗАКЛЮЧЕНИЕ
Проведенные исследования показали, что исключение нуль ядра rot rot -оператора методом деревьев-кодеревьев позволило получить ускорение только в частотном диапазоне 0.00034 - 0.002 Гц по сравнению с использованием конечноэлементной постановки для одного
вектор-потенциала A без исключения нуль ядра rot rot -оператора. Также стоит отметить, что время решения СЛАУ при использовании метода деревьев-кодеревьев оказывается сильно зависящим от направления ребер, включаемых в остовное дерево.
Применение A-V-постановки, в которой для аппроксимации вектор-потенциала A используются реберные базисные вектор-функции, а для скалярного потенциала V - скалярные узловые базисные функции, позволило до 8 раз сократить время решения конечноэлементных СЛАУ, получаемых в задаче трехмерного магнитотеллурического зондирования на частотах от 0.00034 Гц до 0.3 Гц по сравнению с использованием конечноэлементной постановки для
одного вектор-потенциала A .
СПИСОК ЛИТЕРАТУРЫ
[1] Соловейчик Ю.Г., Рояк М.Э., Персова М.Г. Метод конечных элементов для решения скалярных и векторных задач: учеб. пособие. - Новосибирск: Изд-во НГТУ, 2007. - 896 с. («Учебники НГТУ»).
[2] Бердичевский М.Н., Дмитриев В.И. Модели и методы магнитотеллурики. - M.: Научный мир. 2009. -
680 с.
[3] Жданов М.С. Теория обратных задач и регуляризации в геофизике. - М.: Научный мир, 2007. - 712 с.
[4] Спичак В.В. Магнитотеллурические поля в трехмерных моделях геоэлектрики. - М.: Научный мир, 1999. -
204 с.
[5] Домников П.А., Персова М.Г., Соловейчик Ю.Г., Вагин Д.В. Моделирование трехмерных магнитотеллу-рических полей векторным методом конечных элементов и возможности распараллеливания на процессорах с общей памятью // Научный вестник НГТУ. - 2010. - № 3 (40). - С. 87-96.
[6] Персова М.Г., Соловейчик Ю.Г., Домников П.А., Киреева С.В. 3D-моделирование магнитотеллуриче-ских полей с использованием распараллеливания // Актуальные проблемы электронного приборостроения. АПЭП-2010: материалы X международной конференции. - Новосибирск, 22-24 сентября, 2010. - Том 6. - С. 144-149.
[7] Bossavit A. Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements. Academic Press, 1998. - 352 p.
[8] Dular P., Nicolet A., Genon A., Legros W. A discrete sequence associated with mixed finite elements ind its gauge condition for vector potentials // IEEE Transactions on Magnetics. - Vol. 31, no. 3, May 1995. - P. 1356-1359.
[9] Edlinger R., Biro O. A joint vector and scalar potential formulation for driven high-frequency problems using hybrid edge and nodal finite elements // IEEE Trans. Microwave Theory Tech. - Vol. 44. - P. 15-23, Jan. 1996.
[10] Golias N.A., Tsiboukis T.D. Magnetostatics with edge elements: a numerical investigation in the choice of the tree // IEEE Transactions on Magnetics. - Vol. 30, no. 5, September 1994. - P. 2877-2880.
[11] Igarashi H. On the property of the curl-curl matrix in finite element analysis with edge elements // IEEE Transactions on Magnetics. - Vol. 37, no. 5, September 2001. - P. 3129-3132.
[12] Webb J.P. Edge elements and what they can do for you // IEEE Transactions on Magnetics. - Vol. 29, no. 2, March 1993. - P. 1460-1465.
[13] Sogabe T., Zhang S.-L. A COCR method for solving complex symmetric linear systems // Journal of Computational and Applied Mathematics, 199(2007). - P. 297-303.
[14] Zhou L., Walker H. Residual smoothing techniques for iterative methods // SIAM J. Sci. Computing. - Vol. 15, no. 2, 1994. - P. 297-312.
Домников Петр Александрович, аспирант кафедры прикладной математики НГТУ. Основное направление научных исследований - разработка методов решения конечноэлементных систем уравнений при моделировании трехмерных гармонических электромагнитных полей. Имеет 20 публикаций.
E-mail: p_domnikov@mail.ru.
Киреева Светлана Вячеславовна, аспирант кафедры прикладной математики НГТУ. Основное направление научных исследований - разработка вычислительных схем моделирования электромагнитных полей на базе векторного метода конечных элементов для решения задач электромагнетизма. Имеет 5 публикаций.
E-mail: svsib@yandex.ru
Персова Марина Геннадьевна, доктор технических наук, доцент, профессор кафедры прикладной математики НГТУ. Основное направление научных исследований - конечноэлементное моделирование электромагнитных полей в задачах геоэлектрики и электромеханики. Имеет более 80 публикаций, в том числе 2 монографии. E-mail: persova@fpm.ami.nstu.ru.
Соловейчик Юрий Григорьевич, доктор технических наук, профессор, заведующий кафедрой прикладной математики НГТУ. Основное направление научных исследований - конечноэлементное моделирование электромагнитных и тепловых полей. Имеет более 100 публикаций, в том числе 2 монографии.
E-mail: kpmt@fpm.ami.nstu.ru
P.A. Domnikov, S.V. Kireyeva, M.G. Persova, Yu.G. Soloveichik
Finite element modeling of three dimensional magnetotellurics using tree-cotree technique and joint vector and scalar potentials
This paper deals with the three dimensional magnetotellurics modeling using vector finite element method. The comparison of efficiency of computational scheme based on vector potential A using tree-cotree technique and scheme using joint A-V potentials, where vector potential A employs vector edge basis functions and scalar potential V employs scalar nodal basis functions is carried out.
Key words: vector finite element method, tree-cotree method, magnetotelluric sounding.
