CC BY
13Научный вестник НГТУ. -2013. -№ 1(50)
УДК 519.63
Решение трехмерных задач магнитотеллурики в сложных средах с использованием метода конечных элементов
М.Г. ПЕРСОВА, Ю.Г. СОЛОВЕЙЧИК, П.А. ДОМНИКОВ, Т.Г. ШАШКОВА, М.В. АБРАМОВ, Ю.И. КОШКИНА
Рассматриваются вычислительные схемы конечноэлементиого моделирования трехмерных электромагнитных полей в задачах магнитотеллурических зондирований. Приводится описание математического аппарата решения трехмерной задачи на основе векторной, узловой и векторно-узловой математической постановки и примеры численных экспериментов. Оценивается вычислительная эффективность различных методов решения систем линейных алгебраических уравнений.
Ключевые слова: ЗО-моделирование, метод конечных элементов, гармонические электромагнитные поля, магнитотеллурические зондирования.
ВВЕДЕНИЕ
При проведении электроразведочных исследований важнейшую роль играет программно-математический аппарат, используемый для интерпретации полевых данных. При этом получаемые в ходе интерпретации данных геоэлектрические модели могут содержать большое количество трехмерных объектов, и разрабатываемые методы моделирования должны позволять достаточно быстро и точно выполнять расчеты не только для геоэлектрических моделей, содержащих единичные трехмерные объекты, но и для моделей, содержащих несколько десятков трехмерных объектов.
В данной работе будут рассмотрены различные математические постановки и методы решения получаемых в результате конечноэлементной аппроксимации систем линейных алгебраических уравнений (СЛАУ). На примере вычислительных экспериментов, проведенных для геоэлектрических моделей, получаемых на промежуточных этапах интерпретации данных магнитотеллурических зондирований, будет проанализирована вычислительная эффективность различных вычислительных схем.
1. МАТЕМАТИЧЕСКИЙ АППАРАТ МОДЕЛИРОВАНИЯ ТРЕХМЕРНЫХ ЗАДАЧ МАГНИТОТЕЛЛУРИЧЕСКИХ ЗОНДИРОВАНИЙ
Известно, что одной из самых важных составляющих, обеспечивающих вычислительную эффективность конечноэлементных схем, является использование технологии разделения моделируемого поля на нормальную составляющую, создаваемую источником во вмещающей (как правило, горизонтально-слоистой) среде, и на аномальную составляющую, которая описывает поле влияния расположенных во вмещающей среде неоднородностей удельной проводимости (аномальное поле) [1-4].
При этом возможны различные математические постановки, используемые для моделирования гармонического по времени электромагнитного поля, среди которых мы будем рассматривать следующие: постановку для узлового МКЭ [5], постановку для векторного МКЭ [6], а также специальную векторно-узловую постановку.
Статья получена 7 октября 2012 г.
Математическая модель для узлового МКЭ имеет следующий вид [5]
-—ДАа + гюаАа +а8гаау = (а-а" )Еп, (1)
Цо
-а1У(авгааГ) - 1юШу(аАа) = -Ау ((а - а"п ), (2)
в которой под ДА" понимается вектор-функция, определяемая соотношением ДА° =grad(divA0)-rot(rotA0), вектор напряженности аномальной части электрического поля Е° определяется соотношением Е° = -гюА° -gradF, а вектор напряженности нор-
^ п
мальнои части электрического поля Ь определяется из решения соответствующей задачи меньшей размерности в горизонтально слоистой среде. В уравнениях (1), (2) ц(| - магнитная проницаемость вакуума, г - мнимая единица, со - круговая частота поля источника, а -
удельная проводимость трехмерной среды, ст" - удельная проводимость вмещающей (однородной или горизонтально-слоистой) среды.
На внешней границе 5" расчетной области задаются однородные краевые условия на векторный потенциал А" и скалярный потенциал V
= о, = о. (3)
S
Математическая модель для векторного МКЭ имеет вид [6]
—rotrotAа + г®аАа = (а - ап)Еп, (4)
Мо
где аномальная составляющая напряженности электрического поля Е" определяется соотношением Еа = -г® Аа .
И, наконец, рассмотрим математическую модель для векторно-скалярной постановки. Эта модель так же, как и модель (1), (2), является моделью с двумя потенциалами - векторным А" и скалярным V , по в отличие от постановки (1)-(3) для узлового МКЭ в ней вектор-
потенциал А" ищется в виде линейной комбинации векторных базисных функций. Соответствующая математическая модель имеет вид
—rot(rotAа ) + г®а(Аа + gradFa ) = (а-ап )Еп, (5)
Мо
-г® div(agrad Га ) - г^^(аАа ) = -div ((а - ап )Еп ), (6)
и для нее Ёа = -г'ю(Аа + gradF) .
Для решения СЛАУ, получаемых в результате конечноэлементных аппроксимаций рассмотренных выше задач, вообще говоря, могут быть использованы итерационные методы на подпространствах Крылова [7], ориентированные на решение систем с несимметричными разреженными матрицами, например, метод обобщенных минимальных невязок (general minimal residual method- GMRES) [8], метод бисопряженных градиентов (bi-conjugate gradient method - BCG (иногда используется аббревиатура BiCG)) [9], локально-оптимальная схема (ЛОС) [10].
Кроме этих методов в случае, если матрица СЛАУ является комплексно-симметричной (не обязательно эрмитовой), существуют модификации метода сопряженных градиентов (Conjugate Orthogonal Conjugate Gradient Method - COCG) [11] и сопряженных невязок (Conjugate A-Orthogonal Conjugate Residual Method - COCR) [12]. Эти методы имеют формулы, очень похожие на формулы методов сопряженных градиентов и сопряженных невязок для симметричных вещественнозначных матриц.
Для узлового МКЭ вычислительно эффективным является метод, рассмотренный в работе [13] и основанный на минимизации невязки с использованием блочной релаксации (МНБР).
Заметим, что при решении задач магнитотеллурики нормальная составляющая напряженности электрического поля вычислялась из решения одномерной задачи
1 d2
---- А' + таА? = J*.
Mo dz2 * * *
Электрическое поле Е", связанное с вектор-потенциалом А" соотношением Е" = -гсоА", может быть получено в виде
= -т(-А9 sin а, А™ cos а,о),
где а - это угол поворота токов относительно направления с,.
При решении трехмерной задачи в качестве направления Е, обычно используется два ортогональных направления сторонних токов - вдоль осей X и У .
2. РЕЗУЛЬТАТЫ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ
В данной работе приведем результаты вычислительных экспериментов для геоэлектрических моделей, которые были получены на промежуточных этапах интерпретации практических данных.
В качестве первого примера рассмотрим фрагмент геоэлектрической модели, полученной на промежуточном этапе интерпретации данных МТЗ вдоль профиля ЗДВ (Северо-Западный участок) в пределах Дальневосточного Федерального округа.
При использовании узлового МКЭ итерационные методы решения СЛАУ ОМЫЕБ, ВСО, ЛОС и МНБР для данной задачи не сходились.
Проведем сравнение вычислительных затрат, требуемых для моделирования трехмерного электромагнитного поля при использовании А-постановки и А— К-постановки.
Геоэлектрическая модель содержит 63 приповерхностных объекта с удельным сопротивлением от 0.1 Ом-м до 3000 Ом-м. Параметры горизонтально-слоистой среды приведены в табл. 1, на рис. 1 показан ЗБ-вид геоэлектрической модели с 63 объектами.
Вычислительные затраты, потребовавшиеся для решения конечноэлементных СЛАУ, полученных при использовании векторного МКЭ (^4-постановка и А— ^-постановка), приведены в табл. 2.
Данные, приведенные в табл. 2, показывают, что применение А— К-постановки позволило сократить время решения СЛАУ на самой низкой частоте 0.00038 Гц примерно в 19 раз для Х-направления тока и примерно в 40 раз для У-направления тока. На самых высоких частотах 6.25-500 Гц время решения СЛАУ, полученных при использовании А-постановки, было меньшим, чем при использовании А— К-постановки, однако время решения СЛАУ, соответствующих этим частотам, несущественно мало по сравнению со временем решения СЛАУ на низких частотах.
Таблица 1
Параметры вмещающей горизонтально-слоистой среды для задачи с 63 объектами
№ Толщина,м Удельное сопротивление, Ом-м
1 400 1400
2 17000 10000
3 10000 700
4 22000 200
5 60000 300
6 да 150
Рис. 1. ЗБ-вид геоэлектрической модели с 63 объектами
Общее время решения задачи для 21 частоты и двух направлений тока составило 24 ч 0 мин 37 с при использовании А—К-постаповки и 221 ч 22 мин 36 с при использовании ^-постановки.
Проведем вычислительные эксперименты на геоэлектрической модели, полученной на одном из этапов интерпретации данных МТЗ в районе Верхне-Илимпейской параметрической скважины (это одна из обследуемых на наличие минеральных ресурсов площадей Восточной Сибири).
В геоэлектрическую модель включены массивные глубинные аномалии и множество относительно небольших аномалий, расположенных на небольших глубинах (всего 36 объектов). Параметры вмещающей горизонтально-слоистой среды приведены в табл. 3.
Таблш!а 2
Вычислительные затраты для магнитотеллурической задачи с 63 объектами
Ток (Хили У), частота, Гц Векторный МСЭ, СОСЩМЯЗ), А-Г-постановка Векторный МСЭ, СОСЩМКЗ), А-постановка
итерации время итерации время
X, 0.00038 6410 48мин 52 с 293378 15 ч 9 мин 41 с
X, 0.00076 6427 49мин 14 с 247750 12 ч 48мин 12 с
Окончание табл. 2
Ток (Хшти У), частота, Гц Векторный ^Э, COCR(MRS), А-Г-постановка Векторный МСЭ, COCR(MRS), А-постановка
итерации время итерации время
X, 0.00152 5197 39мин 41с 161885 8ч21шн 57 с
X, 0.00305 6206 47мин 18с 124208 6 ч 25мин 8 с
X, 0.0061 6073 46мин 3 с 108480 5 ч 36мин 22 с
X, 0.0122 6445 48мин 50 с 66905 3 ч 27мин 27 с
X, 0.02441 5475 41мин 39 с 45317 2 ч 20мин 31с
X, 0.04882 5713 43мин 22 с 30394 1ч 34мин 14 с
X, 0.19531 4313 32мин 43 с 16414 50мин 53 с
X, 0.39062 3495 26мин 32 с 10697 33мин 10 с
X, 0.78125 2039 15мин 28 с 4881 15мин 8 с
X, 1.5625 1162 8мин 53 с 1858 5мин 45 с
X, 3.125 681 5мин 13 с 784 2мин 25 с
X, 6.25 408 3мин 10 с 566 1мин 45 с
X, 12.5 290 2мин 13 с 311 57 с
X, 50 151 1мин 11с 178 33 с
X, 100 113 54 с 110 20 с
X, 250 102 48 с 109 20 с
X, 500 102 49 с 93 17с
У, 0.00038 7789 59мин 21с 786053 40 ч 37мин 21 с
У, 0.00076 8288 1ч 2 мин 58 с 575295 29 ч 43мин 50 с
У, 0.00152 9563 1ч 12шн 28 с 452159 23 ч 22мин 1 с
У, 0.00305 10544 1ч 19шн 54 с 351491 18 ч 9 мин 53 с
У, 0.0061 11422 1ч 26шн 7 с 267834 13ч50мин29с
У, 0.0122 8944 1ч 7 мин 37 с 206099 10 ч 39мин 3 с
У, 0.02441 9196 1ч 9 мин 28 с 167597 8ч39шн 40 с
У, 0.04882 9406 1ч10шш 46 с 116391 6 ч 0 мин 53 с
У, 0.19531 10395 1ч18шн 52 с 53320 2 ч 45мин 19 с
У, 0.39062 8896 1ч 8 мин 1 с 37143 1ч 55мин 10 с
У, 0.78125 6270 47мин 33 с 32261 1ч 40мин 2 с
У, 1.5625 4816 36мин 25 с 17378 53мин 53 с
У, 3.125 3591 27мин 22 с 10973 34мин 1 с
У, 6.25 1992 15мин 13 с 4569 14мин 10 с
У, 12.5 1107 8мин 27 с 1263 3мин 54 с
У, 50 253 1™ 57 с 299 55 с
У, 100 185 1™ 26 с 197 36 с
У, 250 157 1мин 13 с 173 32 с
У, 500 165 1мин 17с 163 30 с
Таблица 3
Параметры горизонтально-слоистой среды для геоэлектрической моделис 36 объектами
№ Толщина слоя,м Сопротивление слоя, Ом-м
1 60 100
2 490 300
3 100 20
4 650 6
5 2500 500
6 5000 700
7 25000 1000
8 1000 16
9 170000 700
10 50000 500
11 10000 100
12 да 50
ЗБ-вид геоэлектрической модели показан на рис. 2. Конечноэлементная сетка содержит 561759 несогласованных параллелепипедов, 588007 узлов и 1781314 ребер, из них 141399 являются терминальными ребрами и 1639915 - регулярными. Вычислительные затраты, требуемые для решения конечноэлементных СЛАУ, получающихся в результате аппроксимации задачи векторным МКЭ и узловым МКЭ, приведены в табл. 4 и 5 соответственно.
Рис. 2. ЗБ-вид геоэлектрической модели с 36 объектами
Таблица 4
Вычислительные затраты для магнитотеллурической задачи с 36 объектами при использовании векторного МКЭ
Ток (Хшш У), частота, Гц Векторный МКЭ, СОСЩМЯЗ), А—Г-постановка Векторный МКЭ, СОСЩМЯЗ), А-постановка
итерации время итерации время
X, 0.000584 3997 1ч53шш 36 с 222483 50 ч 12мин 40 с
X, 0.001166 4128 1ч 56мин 19 с 174904 39 ч 28мин 23 с
X, 0.002326 4253 2 ч 0 мин 6 с 142183 32 ч 5 мин 19 с
X, 0.004642 4222 1ч 58мин 42 с 108071 24 ч 23мин 24 с
X, 0.009261 4491 2 ч 5 мин 52 с 82885 18ч42мин21 с
X, 0.01848 4579 2 ч 7 мин 29 с 60696 13 ч 41мин 53 с
X, 0.03687 3966 1ч 50мин 34 с 47786 10 ч 47мин 4 с
X, 0.07356 3758 1ч 44мин 40 с 33717 7 ч 36мин 33 с
X, 0.1468 2780 1ч 17шн 45 с 22236 5 ч 1 мин 6 с
X, 0.5843 1900 53мин 2 с 11050 2ч29мин37с
X, 1.166 1375 38мин 34 с 7822 1ч 45мин 55 с
X, 2.326 1198 33мин 27 с 6282 1ч 25мин 3 с
X, 4.642 1021 28мин 40 с 4788 1ч 4 мин 50 с
X, 9.261 1008 28мин 19 с 4070 55мин 6 с
X, 36.87 914 25мин 42 с 2084 28мин 13 с
X, 146.8 614 17мин 21с 958 12мин 58 с
X, 316.2 410 11мин34с 620 8мин 23 с
У, 0.000584 4564 2 ч 8 мин 54 с 237470 53 ч 35мин 36 с
У, 0.001166 4694 2ч12шш 13 с 186039 41ч 59мин 10 с
У, 0.002326 4945 2 ч 18мин47с 142355 32 ч 7 мин 38 с
У, 0.004642 4714 2 ч 12мин 26 с 108758 24 ч 32мин 42 с
У, 0.009261 5151 2ч24мин37с 89414 20 ч 10мин46с
У, 0.01848 5174 2ч24мин51 с 68948 15 ч 33мин 38 с
У, 0.03687 4668 2ч10шн 40 с 50109 11ч 18мин 31с
У, 0.07356 4227 1ч 59мин 7 с 34708 7 ч 49мин 59 с
У, 0.1468 3146 1ч 28мин 6 с 23236 5 ч 14мин 38 с
У, 0.5843 1830 51мин 14 с 12050 2 ч 43мин 10 с
У, 1.166 1381 38мин 45 с 8383 1ч 53мин 30 с
У, 2.326 1242 34мин 42 с 7098 1ч 36шн 6 с
У, 4.642 1047 29мин 19 с 5043 1ч 8 мин 17с
У, 9.261 943 26мин 30 с 4457 1ч 0 гаи 21 с
У, 36.87 821 23мин 7 с 2152 29мин 8 с
У, 146.8 605 17мин 8 с 942 12мин 45 с
У, 316.2 445 12мин 37 с 657 8мин 53 с
Таблица 5
Вычислительные затраты для магнитотеллурической задачи с 36 объектами при использовании узлового МКЭ
Ток (Хшш У), частота, Гц Узловой МКЭ,^НБР Узловой МСЭ, 0МКЕ8(50)
итерации время итерации время
X, 0.000584 4 5мин 58 с 239 5 ч 5 мин 39 с
X, 0.001166 6 10мин 8 с 276 5ч52шш 58 с
X, 0.002326 8 12мин 42 с 328 6 ч 59мин 28 с
X, 0.004642 14 22мин 8 с 493 10ч30мин29с
X, 0.009261 29 45мин 25 с >1000 >21ч
У, 0.000584 4 6мин 31с 262 5 ч 35мин 3 с
У, 0.001166 6 9мин 44 с 286 6 ч 5 мин 45 с
У, 0.002326 9 14мин 24 с 357 7 ч 36мин 33 с
У, 0.004642 15 24мин 3 с 712 15 ч 10мин 33 с
У, 0.009261 29 46мин 17с >1000 >21ч
При использовании узлового МКЭ и метода МНБР решение СЛАУ удалось получить только для нижних 5 частот от 0.000584 Гц до 0.009261 Гц. Метод ОМЫЕБ позволил получить решение только для четырех нижних частот. Для остальных частот метод ОМЫЕБ не смог понизить невязку до требуемого уровня за 1000 итераций (что соответствовало примерно 21 ч времени).
Общее время решения задачи (на всех частотах), полученное для А -постановки и Л—К-постановка, а также для их комбинации с методом МНБР, приведено в табл. 6 (на 1 ядре и 10 ядрах). Распараллеливание выполнялось по гармоникам.
Таблица 6
Время решения задачи на всех частотах
Метод решения Время (1 ядро) Время (10 ядер)
Векторный МЕСЭ, СОСЯ, 4-постановка 439 ч 16мин 61ч 44мин
Векторный МКЭ, СОСЯ, 4-Г-постановка 46 ч 20мин 4 ч 45мин
Комбинирование методов ^НБР и СОСЯ (4-постановка) 105 ч 16мин 17 ч 34мин
Комбинирование методов ^НБР и СОСЯ (4-Т-постановка) 28 ч 26мин 3 ч 37мин
ЗАКЛЮЧЕНИЕ
Результаты вычислительных экспериментов показали, что время решения задач магнитотеллурики в средах с большим числом ЗБ-объектов на низких и средних частотах при использовании А— ^-постановки на порядок меньше, чем при использовании А -постановки. Наиболее выгодным по вычислительным затратам оказалось комбинирование узлового МКЭ с решателем МНБР на низких частотах и векторного МКЭ в Л-К-постаповке с решателем СОСЯ на средних и высоких частотах.
Выявленные в результате проведенных исследований наиболее эффективные постановки и методы позволяют достаточно устойчиво и за приемлемое время получать решение самых сложных трехмерных задач магнитотеллурики и могут быть использованы для интерпретации данных магнитотеллурических зондирований.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг.
СПИСОК ЛИТЕРАТУРЫ
[1] ПерсоваМ.Г. Компьютерное моделирование геоэлектромагнитных полей в трехмерных средах методом конечных элементов /М.Г. Персова, Ю.Г. Соловейчик, Г.М. Тригубович// Физика Земли, 2011. -№ 2. - С. 3-14.
[2] Соловейчик Ю.Г. Метод конечных элементов для решения скалярных и векторных задач / Ю.Г. Соловейчик, М.Э. Рояк, М.Г. Персова. - Новосибирск: НГТУ, 2007. - 896 с.
[3] Соловейчик Ю.Г. Моделирование нестационарных электромагнитных полей в трехмерных средах методом конечных элементов / Ю.Г. Соловейчик, М.Э. Рояк, B.C. Моисеев, Г.М. Тригубович // Физика Земли, 1998. - № 10. -С. 78-84.
[4] Badea Е.А. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials / E.A. Badea, M.E. Everett, G.A. Newman, O. Biro // Geophysics, 2001. - Vol. 66. - № 3. - P. 786-799.
[5] Рояк М.Э. Конечноэлементное моделирование трехмерных гармонических электромагнитных полей в задачах аэроэлектроразведки кимберлитовых трубок / М.Э. Рояк, С.Х. Рож, Ю.Г. Соловейчик, Г.М. Тригубович // Сибирский журнал индустриальной математики. - 1998. - Т. 1, № 2. - С. 154-168.
[6] Соловейчик Ю.Г. Метод конечных элементов для решения скалярных и векторных задач: учеб. пособие / Ю.Г. Соловейчик, М.Э. Рояк, М.Г. Персова. - Новосибирск: Изд-во НГТУ, 2007. - 896 с.
[7] Ильин В.П. Методы и технологии конечных элементов / В.П. Ильин. - Новосибирск: Изд-во ИВМиМГ, 2007.-371 с.
[8] SaadY. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems / Y. Saad, M. Schultz // SIAM J. Sci. Comput. - 1986- Vol. 7. - P. 856-869.
[9] Saad Y. Iterative methods for sparse linear systems - SIAM, Philadelphia, - 2003. - 2nd ed. - 528 Pp.
[10] Soloveichik Y.G. Iterative method for solving finite element systems of algebraic equations / Y.G. Soloveichik // Computers & Mathematics with Applications. - Vol. 33, Issue 6. - 1997. - P. 87-90.
[11] Van der Vorst Н.А. A Petrov-Galerkin type method for solving Ax=b, where A is symmetric complex / H.A. Vorst van der, J.B.M. Melissen // IEEE Transaction on Magnetics. - Vol. 26. - № 2. - 1990. - P. 706-708.
[12] Sogabe T. A COCR method for solving complex symmetric linear systems / T. Sogabe, S.-L. Zhang // Journal of Computational and Applied Mathematics. - № 199. - 2007. - P. 297-303.
[13] Домников П.А. Метод решения систем уравнений, возникающих при конечноэлементной аппроксимации гармонических по времени электромагнитных полей / П.А. Домников // Сборник научных трудов НГТУ. - 2009. -№2(56).-С. 41-46.
Персова Марина Геннадьевна, доктор технических наук, профессор, профессор кафедры прикладной математики НГТУ. Основное направление научных исследований: конечноэлементное моделирование электромагнитных полей в задачах геоэлектрики и электромеханики. Имеет более 80 публикаций, в том числе 2 монографии. E-mail [email protected]
Соловейчик Юрий Григорьевич, доктор технических наук, профессор, заведующий кафедрой прикладной математики НГТУ. Основное направление научных исследований: конечноэлементное моделирование электромагнитных и тепловых полей. Имеет более 100 публикаций, в том числе 2 монографии. E-mail [email protected]
Домников Петр Александрович, аспирант кафедры прикладной математики НГТУ. Основное направление научных исследований: разработка методов решения СЛАУ, полученных в результате конечноэлементной аппроксимации в задачах электромагнетизма. Имеет 34 публикации. E-mail [email protected]
ГПашкова Татьяна Геннадьевна, аспирант кафедры прикладной математики НГТУ. Основное направление научных исследований: конечноэлементное моделирование электромагнитных полей в задачах геоэлектрики. Имеет 6 публикаций. E-mail [email protected]
Абрамов Мшат Владимирович, кандидат технических наук, доцент кафедры прикладной математики НГТУ. Основное направление научных исследований: конечноэлементное моделирование электромагнитных полей в задачах геоэлектрики. Имеет 20 публикаций. E-mail [email protected]
Кошкина Юлия Игоревна, магистрант кафедры прикладной математики НГТУ. Основное направление научных исследований: решение обратных задач геоэлектроразведки. Имеет 1 публикацию. E-mail [email protected]
Persova M.G., Soloveichik Yu.G., Domnikov P.A., Shashkova T.G.,Abramov M.V., Koshkina Yu. I.
Solution of three-dimensional magnetotelluricsproblems using thefinite element method in complicated media
Computational schemes for finite element modelling of 3-D electromagnetic fields in magnetotelluric sounding problems are considered. The description of mathematical apparatus for solving the 3-D problem on the basis of the vector, node and vector-node mathematical statement, and examples of numerical experiments are presented. The computational efficiency ofvarious methods used for solving the systems oflinear algebraic equations is estimated.
Key words: 3-D modeling, finite element method, harmonic electromagnetic fields, magnetotelluric soundings.
