автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Построение дискретных аналогов и разработка методов решения конечноэлементных систем уравнений при моделировании гармонических электромагнитных полей

На правах рукописи

Домников Петр Александрович

ПОСТРОЕНИЕ ДИСКРЕТНЫХ АНАЛОГОВ И РАЗРАБОТКА МЕТОДОВ РЕШЕНИЯ КОНЕЧНОЭЛЕМЕНТНЫХ СИСТЕМ УРАВНЕНИЙ ПРИ МОДЕЛИРОВАНИИ ГАРМОНИЧЕСКИХ ЭЛЕКТРОМАГНИТНЫХ ПОЛЕЙ

Специальность 05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата технических наук

2 9 НОЯ /1Л2

Новосибирск - 2012

005055932

005055932

Работа выполнена в Федеральном государственном бюджетном образовательном учреждении высшего профессионального образования «Новосибирский государственный технический университет»

Научный руководитель: доктор технических наук, профессор

Соловейчик Юрий Григорьевич

Официальные оппоненты: Селезнев Вадим Александрович, доктор фи-

зико-математических наук, профессор, Новосибирский государственный технический университет, заведующий кафедрой высшей математики

Трусов Алексей Андреевич, кандидат технических наук, главный геофизик ЗАО «ГНПП Аэрогеофизика»

Ведущая организация: Федеральное государственное бюджетное

учреждение науки Институт вычислительной математики и математической геофизики СО РАН

Защита состоится « 13 » декабря 2012 г. в 16-00 часов на заседании диссертационного совета Д 212.173.06 при Федеральном государственном бюджетном образовательном учреждении высшего профессионального образования «Новосибирский государственный технический университет» по адресу: 630092, г. Новосибирск, пр. К. Маркса, 20.

С диссертацией можно ознакомиться в библиотеке Новосибирского государственного технического университета.

Автореферат разослан «11» ноября 2012 года.

Ученый секретарь

диссертационного совета ¿Щъ-- Чубич Владимир Михайлович

Общая характеристика работы

Задачи моделирования трехмерных геоэлектромагнитных полей в неоднородных средах относятся к одним из наиболее вычислительно трудоемких. Уменьшение вычислительных затрат и повышение точности решения этих задач, несомненно, является актуальной задачей, поскольку от этого напрямую зависит качество, а часто и сама возможность адекватного восстановления геоэлектрического строения среды.

Основными методами, применяемыми для численного моделирования трехмерных геоэлектромагнитных полей, являются метод интегральных уравнений (МИУ), метод конечных разностей (МКР) и метод конечных элементов (МКЭ). При этом известно, что МИУ, в разработку которого существенный вклад внесли В.И. Дмитриев, М.С. Жданов, Z. Xiong, Р.Е. Wannamaker и др., является вычислительно эффективным при расчетах полей в средах с небольшим количеством трехмерных объектов, с ростом же количества трехмерных неоднородностей (что более соответствует реальным условиям) этот метод начинает существенно уступать по вычислительной эффективности таким сеточным методам, как МКР (для задач геоэлектрики он разрабатывался В.Л. Дру-скиным, JI.A. Книжнерманом, D.L. Alumbaugh, R.L. Mackie, G.A. Newman, T. Wang и др.) и МКЭ, в развитие которого для решения задач геоэлектромагнетизма значимый вклад внесли М.Г. Персова, Ю.Г. Соловейчик, L. Changsheng, М.Е. Everett, Y. Mitsuhata, M.J. Nam, X.Y. Sun и др. МКР и МКЭ основаны на аппроксимациях дифференциальных уравнений в частных производных. При этом по вычислительной эффективности МКР, в свою очередь, уступает МКЭ из-за необходимости использования в МКР прямоугольных регулярных сеток -для получения решения требуемой точности это приводит к дискретным аналогам гораздо большей размерности.

В принципе, МКР и МКЭ уже достаточно давно применяются для решения трехмерных задач электроразведки. Основным преимуществом этих методов перед МИУ является разреженность матриц систем уравнений, получаемых в результате аппроксимаций соответствующих трехмерных задач. Но вместе с тем МКР и МКЭ требуют включения в расчетную область помимо самих трехмерных объектов довольно большого пространства вокруг них, и при этом для достижения хорошей точности необходимы достаточно подробные сетки, причем мелкие ячейки нужны не только в ЗБ-объектах, но и в окружающем их пространстве, особенно вблизи источников поля и приемников. В результате, при решении задач с контролируемыми источниками, если ЗБ-объекты в них дают относительно слабые отклики, МКР и МКЭ требуют очень высоких вычислительных затрат на получение решения с приемлемой точностью.

Однако разрабатываемые в последние годы модификации МКЭ с включением в него возможности выделения поля простой структуры (что является стандартным для МИУ и не использовалось до самого последнего времени в МКЭ) - так называемого нормального поля или поля вмещающей среды — существенно расширили его возможности при решении трехмерных задач геоэлектрики. Присущие МКЭ возможности учета геометрии расчетной области любого уровня сложности и его гибкость при построении сеток в сочетании с

технологиями, использующими выделение нормального поля, позволяют сделать этот метод мощным и удобным инструментом при решении трехмерных задач геоэлектромагнетизма.

При использовании МКЭ для моделирования трехмерных геоэлектромагнитных полей основная вычислительная нагрузка, как правило, приходится на решение систем линейных алгебраических уравнений (СЛАУ), получаемых в результате конечноэлементной аппроксимации исходной краевой задачи.

При решении трехмерных задач, особенно векторных, конечноэлемент-ные аппроксимации приводят к СЛАУ с разреженными матрицами очень большого размера с довольно большой шириной ленты, поэтому для решения таких СЛАУ применяют, в основном, итерационные методы, которые рассмотрены, например, в работах О. Axelsson, C.G. Broyden, Y. Saad, H.Van der Vorst, В.П. Ильина и др. При этом скорость сходимости различных итерационных методов может сильно (в десятки, а иногда и в сотни раз) различаться, причем при решении конкретной практической задачи, как правило, заранее неизвестно, какой из методов решения может оказаться наиболее быстрым. Поэтому исследование различных методов решения СЛАУ, получаемых в результате конечно-элементной аппроксимации трехмерных задач, также является актуальной задачей.

При моделировании трехмерных векторных полей в методе конечных элементов наиболее широко распространенными являются следующие два подхода. Первый из них (узловой или скалярный МКЭ) состоит в том, что каждая из компонент искомого поля ищется в виде линейной комбинации скалярных базисных функций (которые обычно связывают с узлами конечноэлементной сетки). Второй подход (векторный МКЭ) позволяет строить математические модели непосредственно в терминах векторных переменных (используются векторные базисные функции). В условиях, когда токами смещения можно пренебречь (моделирование электромагнитных полей на не очень высоких частотах), можно использовать постановки как для векторных, так и для скалярных базисных функций, и в зависимости от решаемой задачи преимущество может иметь как одна, так и другая постановка.

Отметим, что даже при использовании технологии выделения поля вычислительные затраты, требуемые для расчета трехмерных полей в сложных средах, содержащих много трехмерных включений, могут оставаться довольно большими. Дополнительный эффект сокращения вычислительных затрат может быть получен за счет оптимизации сеток с обоснованным изменением шагов и устранением так называемых "лишних" узлов из регулярных сеток. В этом случае в сетках могут появляться терминальные узлы, и необходима разработка специальных подходов к построению дискретных аналогов. Кроме того, вычислительные затраты при решении задач геоэлектромагнетизма с использованием МКЭ довольно сильно зависят от размера расчетной области. Ее границы должны быть удалены от мест регистрации сигнала настолько, чтобы краевые условия практически не влияли на получаемые в приемниках значения, поэтому важной является задача определения удаленных границ расчетной области (границ так называемого большого объема или "бака"). Реализация возможно-

ста определения минимально допустимого размера "бака" позволяет дополнительно сократить размер соответствующих дискретных аналогов и, соответственно, уменьшить вычислительные затраты.

В целом, все рассмотренные проблемы, связанные со снижением вычислительных затрат при конечноэлементном моделировании трехмерных геоэлектромагнитных полей, решению которых посвящена данная диссертационная работа, являются отдельными аспектами основной научной проблемы создания эффективных методов численного моделирования трехмерных гармонических по времени электромагнитных полей.

Цели и задачи исследования

1. Построение эффективных конечноэлементных аппроксимаций для решения задачи моделирования трехмерных гармонических по времени электромагнитных полей. Для достижения данной цели сформулированы и решены следующие задачи:

- разработка и программная реализация конечноэлементных вычислительных схем моделирования гармонических по времени электромагнитных полей на трехмерных несогласованных сетках на основе узлового и векторного МКЭ с использованием выделения основной части поля;

- построение дискретных аналогов при использовании метода деревьев-кодеревьев и совместного использования скалярного и векторного потенциалов в векторном МКЭ на трехмерных сетках; реализация разработанных вычислительных схем в программном комплексе.

2. Разработка и реализация эффективных процедур решения конечноэлементных СЛАУ, возникающих при конечноэлементной аппроксимации трехмерных гармонических по времени электромагнитных полей. Для достижения данной цели сформулированы и решены следующие задачи:

- экспериментальное исследование существующих методов решения СЛАУ, возникающих при использовании различных конечноэлементных постановок для моделирования гармонических электромагнитных полей;

- разработка методов решения СЛАУ, получаемых при использовании узлового МКЭ для моделирования трехмерных гармонических электромагнитных полей;

- программная реализация методов решения конечноэлементных СЛАУ, получаемых в результате аппроксимации трехмерных гармонических по времени электромагнитных полей; сравнение производительности разработанных программ с существующими математическими библиотеками.

На защиту выносятся

1. Математические постановки и схемы конечноэлементного моделирования трехмерных гармонических по времени электромагнитных полей с использованием узлового и векторного МКЭ с выделением части поля на сетках с терминальными узлами.

2. Методы решения СЛАУ, возникающих при моделировании трехмерных электромагнитных полей на высоких частотах с использованием векторного МКЭ и на низких частотах с использованием узлового МКЭ.

3 Программные модули конечноэлементного моделирования трехмерных

' геоэлектромагнитных полей и решения СЛАУ в подсистеме ЭР-МТЗ программного комплекса ОеоЕМ.

4. Результаты использования разработанных методов и программного комплекса для решения задач геоэлектромагнетизма. Достоверность результатов

Вычислительные схемы конечноэлементного моделирования трехмерных электромагнитных полей протестированы путем сравнения с решением задач меньшей размерности, с решениями, полученными другими авторами, а также путем сравнения решений одних и тех же задач с использованием представлений электромагнитного поля через различные потенциалы в векторном и в узловом мкэ.

Научная новизна

1. Разработана и реализована вычислительная схема для моделирования трехмерных магнитотеллурических полей в математической постановке для векторного МКЭ с использованием векторного и скалярного потенциалов, технологии выделения поля вмещающей горизонтально-слоистой среды и нерегулярных сеток с терминальными узлами.

2. Разработана и реализована вычислительная схема для моделирования трехмерных гармонических геоэлектромагнитных полей с использованием метода деревьев-кодеревьев и технологии выделения части поля.

3. Разработан новый метод решения СЛАУ, возникающих при моделировании трехмерных гармонических по времени электромагнитных полей на низких частотах с использованием узлового МКЭ - метод минимизации невязки с использованием блочной релаксации (МНБР). Практическая значимость

Предлагаемые в данной работе способы построения дискретных аналогов и методы решения СЛАУ реализованы в подсистемах ЭР-МТЗ и ЭРИАН программного комплекса СеоЕМ, который использовался при проектировании поисково-разведочных работ и интерпретации практических данных. Публикации

По результатам выполненных исследований опубликовано 29 работ, в том числе 7 статей в журналах, входящих в перечень изданий, рекомендуемых Высшей аттестационной комиссией, 2 статьи в сборнике научных трудов, 18 публикаций в материалах конференций, 2 свидетельства о государственной регистрации программы для ЭВМ. Апробация работы

Основные результаты работы были представлены и докладывались на: международной научно-практической конференции "Геобайкал 2012" (Иркутск, 2012); всероссийской школе-семинаре им. М.Н. Бердичевского и Л.Л. Ваньяна по электромагнитным зондированиям Земли, ЭМЗ-2011 (Санкт-Петербург, 2011); рабочем семинаре "Наукоемкое программное обеспечение" в рамках Ершовской конференции по информатике (РБШ) (Новосибирск, Академгородок, 2011); IX и X международной конференции "Актуальные проблемы электронного приборостроения" АПЭП-2008, АПЭП-2010 (Новосибирск, 2008,

2010); всероссийской конференции по вычислительной математике КВМ-2011 (Новосибирск, Академгородок, 2011); научной конференции Лаборатории НГУ-Интел "Технологии высокопроизводительных вычислений" (Новосибирск, Академгородок, 2009, 2010); IV международной научно-технической конференции молодых специалистов, аспирантов и студентов "Математическое и компьютерное моделирование естественнонаучных и социальных проблем" (Пенза, 2010); всероссийской научной конференции молодых ученых "Наука. Технологии. Инновации" (Новосибирск, 2008, 2009, 2010); российской научно-технической конференции "Информатика и проблемы телекоммуникаций" (Новосибирск, 2007, 2008, 2010, 2011); VI Всероссийской научно-практической конференции студентов, аспирантов и молодых ученых "Технологии Microsoft в теории и практике программирования" (Томск, 2009).

Личный вклад

Разработан метод решения СЛАУ, возникающих при моделировании гармонических по времени трехмерных электромагнитных полей на низких частотах с использованием узлового МКЭ. Разработаны, реализованы и апробированы программные модули конечноэлементного моделирования и решения СЛАУ в программном комплексе GeoEM. Разработан и реализован алгоритм подбора размера расчетной области для решения трехмерных задач магнитотеллурики. Выполнена компьютерная реализация методов решения конечноэлементных СЛАУ в многоядерных системах с общей памятью и реализация запуска задач магнитотеллурических зондирований по частотам в системах с распределенной памятью. Разработаны и исследованы методы построения дискретных аналогов при моделировании трехмерных электромагнитных полей с использованием метода деревьев-кодеревьев.

В совместных публикациях автору принадлежат следующие результаты. В работах [4, 26, 27] автором проведены расчеты трехмерных электромагнитных полей, в работах [2, 5, 16, 23, 28, 29] автором разработаны методы решения конечноэлементных СЛАУ, в работах [6, 25] автором разработаны и реализованы модули распараллеливания расчетов трехмерных электромагнитных полей и выполнены вычислительные эксперименты, в [1, 15, 24] проведены сравнения итерационных методов решения СЛАУ.

Работа выполнялась при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы (автор принимал участие в 9 проектах, в том числе в 6 проектах являлся руководителем) и ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2013 годы» (в 3 проектах). Работа выполнялась в рамках государственных заданий высшим учебным заведениям в части проведения НИР (шифр заявки 8.874.2011), а также была поддержана грантом Фонда содействия отечественной науке в рамках программы «Лучшие аспиранты РАН» в 2010 г.

Структура работы

Диссертационная работа состоит из введения, четырех глав, заключения, списка использованных источников (133 наименования), приложения и содержит 40 рисунков и 52 таблицы. Общий объем работы -191 страница.

Краткое содержание работы

В первой главе проведен обзор научной литературы по теме диссертационного исследования и приведены математические модели в виде дифференциальных уравнений и систем дифференциальных уравнений, использующиеся для моделирования трехмерных гармонических по времени геоэлектромагнитных полей и ориентированные на использование векторного и узлового МКЭ.

Основными методами, применяемыми для численного моделирования трехмерных геоэлектромагнитных полей, являются метод интегральных уравнений (МИУ), метод конечных разностей (МКР) и метод конечных элементов (МКЭ). Первые значимые результаты при расчетах трехмерных геоэлектромагнитных полей были получены с помощью МИУ. Однако МИУ [Д.Б. Авдеев, В.И. Дмитриев, М.С. Жданов, A. Abubakar, Е. Liu, Z. Xiong, P.E. Wannamaker, В. Singer] с ростом аномальной области (а это характерно для большинства практических задач) резко теряет свои преимущества. Это связано с тем, что аппроксимация задачи по методу интегральных уравнений приводит к системе линейных алгебраических уравнений (СЛАУ) с плотной матрицей, и МИУ начинает существенно уступать по вычислительной эффективности таким сеточным методам, как МКР [Ив.М. Варенцов, B.JI. Друскин, JI.A. Книжнерман, D.L. Alumbaugh, R.L. Mackie, G.A. Newman, T. Wang и др.] и МКЭ [М.Г. Персова, Ю.Г. Соловейчик, L. Changsheng, М.Е. Everett, Y. Mitsuhata, M.J. Nam, X.Y. Sun и др.], основным преимуществом которых является разреженность матриц систем уравнений, получаемых в результате аппроксимаций соответствующих трехмерных задач.

В научных публикациях довольно часто обсуждаются проблемы, связанные с медленной скоростью сходимости итерационных методов и с предобу-словливанием при решении задач моделирования электромагнитных полей в технических устройствах с использованием векторного МКЭ (ВМКЭ) [A. Bossavit, P. Dular, О. Biro, N. Golias, H. Igarashi и др.]. При этом отмечается, что одной из возможных причин медленной сходимости итерационных методов при решении СЛАУ, получаемых в результате ВМКЭ-аппроксимации задач электромагнетизма, является наличие большого нуль-ядра у оператора rot rot (состоящего из функций вида gradç?). Для устранения нуль-ядра rotrot-оператора существует технология деревьев-кодеревьев (tree-cotree technique) [A. Bossavit, P. Dular, R. Dyczij-Edlinger и др.]. В научной литературе наблюдаются серьезные разногласия между разными авторами по поводу эффективности использования метода деревьев-кодеревьев. В одной части работ приводятся результаты успешного применения метода tree-cotree для решения практических задач, и в то же время Н. Igarashi сделано наблюдение, что применение метода tree-cotree может ухудшать спектральные свойства матрицы СЛАУ и, таким образом, применение метода tree-cotree может и не приводить к уменьшению времени решения СЛАУ. Применение метода tree-cotree для задач геоэлектромагнетизма, где матрицы СЛАУ обладают несколько другими свойствами, в научной литературе практически не рассмотрены.

При всем многообразии используемых для численного моделирования трехмерных геоэлектромагнитных полей подходов в них практически не рассматриваются задачи с большим числом ЗО-объектов, что требуется для обработки реальных данных геоэлектрических зондирований. Кроме того, в подходах, основанных на МКР и МКЭ, как правило, не используется технология выделения поля горизонтально-слоистой среды, которая во многих случаях позволяет уменьшить время решения трехмерных задач более чем на порядок.

В работах Ю.Г. Соловейчика в середине 1990-х годов была разработана модификация МКЭ с включением в него возможности выделения поля простой структуры (что является стандартным для МИУ и не использовалось в МКЭ) — так называемого нормального поля или поля вмещающей среды, что позволило без снижения точности численного решения использовать достаточно грубые дискретизации по пространству и, соответственно, существенно расширило возможности МКЭ при решении трехмерных задач геоэлектрики. Позднее эта схема (для контролируемых источников) была независимо предложена в работах зарубежных исследователей [Е.Badea, M.Everett (2001 г.); S. Mukheijee, М. Everett (2011 г.)].

Приведем математические модели для описания трехмерных гармонических электромагнитных полей, основанные на использовании технологии выделения поля горизонтально-слоистой среды. Вектор-потенциал магнитного поля А (являющийся комплекснозначной величиной) представляется в виде суммы А = Лп +Аа, где вектор-функция А" (нормальное поле) описывает поле в горизонтально-слоистой среде, а вектор-функция Аа представляет собой добавочное поле от трехмерных неоднородностей. Добавочное поле, описываемое вектор-потенциалом Аа, может быть найдено из решения векторного уравнения

—rotrot Аа +icoaAa =(сг-ап)Ёп (1)

МО

с однородными краевыми условиями первого рода (для касательных составляющих А), заданными на границе расчетной области. В (1) a = a(x,y,z) — функция, описывающая удельную проводимость среды с включенными в нее

трехмерными неоднородностями, Ё" - -icoAn - нормальная составляющая напряженности электрического поля во вмещающей горизонтально-слоистой среде (пересчитывается в трехмерную сетку). Индукция магнитного поля Ва и напряженность электрического поля Еа, соответствующие аномальному полю, определяются соотношениями: Ва = rot Аа, Ёа =—ia>Aa.

Уравнение (1) может быть решено с использованием векторного метода конечных элементов (ВМКЭ). Математическая модель электромагнитного поля (1) при использовании ВМКЭ позволяет решать самые сложные задачи геоэлектромагнетизма. Она корректно описывает геоэлектромагнитное поле в ситуациях, когда среда содержит неоднородности с измененными электрическими и магнитными свойствами. На основе этой модели можно моделировать элек-

тромагнитные поля как в непроводящих, так и в проводящих и в смешанных (т.е. содержащих непроводящие и проводящие подобласти) средах. Рассмотрим две вариационные постановки:

— \то1Аа -xotipán + ia\<jAa \{а-а")Ё"-y/áQ., (2)

п Q

которую далее будем называть А -постановкой, и

— f rot Аа - rotvj/dQ + ico¡ a(Aa + gradF")• ydí2 = J (o -о")Ё" ■ y dft, (3) n o £1

ico J a(Aa + grad Va) • grad <pd£l = J (o - а" )Ё" ■ grad cpáQ., (4)

a n

которую далее будем называть A-V-постановкой. В (2)-(4) вещественная и мнимая части решения Аа и вещественная пробная вектор-функция у/ принадлежат пространству Я0(rot.fi), а скалярная функция Va - пространству Я1 функций с суммируемыми с квадратом производными. Напряженность аномального электрического поля и индукция аномального магнитного поля для постановки (3)-(4) вычисляются по формулам Ёа =-ico(Aa + gradK"), Ba = T0tAa.

В случае узлового МКЭ система дифференциальных уравнений для вычисления трехмерного поля (влияния ЗО-неоднородностей) будет иметь вид

__1_да° + icocrAa + agradé =(сг- а")Ё", ■ "0 (5)

-írtv(crgradFa) - icodiv(o-^a) = -div[{a - <т")£"|,

где A = graddiv-rotrot - оператор Лапласа. В данном случае Ba=rotAa, Ёа =-icoАа - gradV". При этом потенциалы Аа и Va аппроксимируются узловыми базисными функциями (для вектор-потенциала Аа имеются ввиду три его компоненты).

Рассмотренные математические модели для вычисления трехмерной части поля (после выделения поля источника во вмещающей горизонтально-слоистой среде) являются универсальными для любого гармонического по времени типа источника поля, который определяет поле Ё" в правых частях соответствующих уравнений. Они могут применяться для моделирования трехмерных магнитотеллурических полей, моделирования трехмерных геоэлектромагнитных полей в гармоническом режиме в задачах аэроэлектроразведки с петлевым источником и в задачах геоэлектроразведки (наземной и на шельфе) с использованием источника в виде ГЭЛ.

Во второй главе описано построение дискретных аналогов для вариационных постановок, представленных в первой главе, на основе узлового и век-

торного МКЭ и рассмотрены методы решения соответствующих конечноэле-ментных СЛАУ.

В данной работе использовалось разбиение расчетной области Q на конечные элементы в форме параллелепипедов и шестигранников. В случае узлового МКЭ использовались трилинейные базисные функции из пространства Я(П), связанные с узлами конечноэлементной сетки.

Матрица конечноэлементной СЛАУ Aq= f, полученная при использовании узлового МКЭ, обладает разреженной блочной структурой. Для ее решения этой СЛАУ использовались методы на подпространствах Крылова, предназначенные для решения систем с несимметричными матрицами: метод GMRES [Y. Saad], метод бисопряженных градиентов (Bi-conjugate gradient method-BCG) [H. Van der Vorst], локально-оптимальная схема (ЛОС) [Ю.Г. Соловейчик]. Также был разработан специальный метод, названный методом минимизации невязки с использованием блочной релаксации (МНБР). Алгоритм метода МНБР выглядит следующим образом. Выбрать начальное приближение

, Re Л, Re A, Re A? Re К ImAx ImAy Im Л, 1тКчГ /дч

4(0) =(4(0) >4(0) >4(0) >4(0) >4(0) >1(0) >4(0) >4(0) ) • (°)

Вычислить r0= f - Aq§- (7)

Для к = 0,1,...,£тах (8) если к - четное:

<&) = С-'(/Кс К + + Df^y + ОЦ^ ); (9)

4^2) -Dzql&,2)+Cql™A*y, (10)

4% = S-V/1^' - + Сс/^У ); (11) 4^1f2) = B-\fReA*(12)

ImAx ImAx ЬnAy 1тЛу imJ ЬлА, ImV ImV „,..

4(M/2) = 4(к) > 4{k+ih) = 4(к) ' 4(k+lh) = 4(k) • 4(k+l/2) = 4(k) i (13) если к — нечетное:

4'&)-G-\flmV-D^ -Df^y-DTzq^y, (14)

4«¿h) = B-\flmA:-Dzq^m)-Cq^y, (15)

Im Ay j Im Ay тщр- ,л/-х

4(к+\/2) (/ -°У4(к+\12)-С4{к) У> (16)

= -DxclT^-C^*)- (17)

Re Л* Re/lv Re/! у ReAy яел, RеА, Re К Re У

4(к+1/2)=4(к) >4(k+U2)=4(k) >4(k+\h)=4(k) > 4(к+\П) =4(к) > (1»)

г(к)=Ч(к+1/2) - <¡(k); (!9)

„ .. (г(к)'Аг(к)) .

ак -»

(А7(к)'Аг(к))

если а/с <£[ — выход; Ч(к+\)=Ч(к)+ак2(к)>

если ||г(А+1) | / |г(0)|| <е2- выход, где В,с,0х,0у,02,0 - вещественные разреженные матрицы, структурными элементами матрицы А:

(21) (22)

(23)

(24)

являющиеся

А =

' в о о Dx -С о о о

о в о Dy О -с о о

о о в DZ О о -с о

о о о G -Dl -DÇ -Dl О

с о о О В о О Dx

о с о О О в о Dy

о о с О о о в DZ

к $ Dl О о о О G

(25)

В случае векторного МКЭ использовались базисные вектор-функции из пространства H(rot,Cï), связанные с ребрами конечноэлементной сетки. Для векторного МКЭ рассматривается построение дискретных аналогов для А-постановки (2), А - F-постановки (3)-(4) и с применением метода деревьев-кодеревьев. При использовании метода деревьев-кодеревьев был рассмотрен случай построения дерева во всей расчетной области и только в непроводящей подобласти. Для каждого из этих двух случаев рассматриваются варианты построения остовного дерева, в которых внутри области в дерево включались ребра, направленные только вдоль одной из координатных осей: либо X (X-дерево), либо Y (K-дерево), либо Z (Z-дерево), а также дерево с произвольно ориентированными ребрами (ATZ-дерево).

Для решения СЛАУ, полученных при использовании векторного МКЭ, использовались методы GMRES, BCG, ЛОС и модификации метода сопряженных градиентов - Conjugate Orthogonal Conjugate Gradient Method (COCG) [H. Van der Vorst] и сопряженных невязок — Conjugate A-Orthogonal Conjugate Residual Method (COCR) [T. Sogabe], предназначенные для решения СЛАУ с комплексно-симметричными матрицами.

Для решения практических задач особенно эффективным является использование несогласованных сеток с локальными сгущениями вблизи объектов. В данной работе была использована технология согласования для гексаэд-ральных конечных элементов с помощью матрицы перехода. Суть ее состоит в

следующем. Пусть согласованные базисные функции <pf, i = \...nc, выражаются

через базисные функции <р"с,i = l...nnc, построенные стандартным способом как сшивка локальных базисных функций конечных элементов, с помощью соотношения фс = Т<р"с, где Т - прямоугольная пс х и матрица перехода от базиса {<р"с} к базису {tpf}. Тогда матрица конечноэлементной СЛАУ преобразуется к

форме Ас = ТА"СТТ, где Апс - матрица, собранная из локальных матриц несогласованных конечных элементов стандартным способом. Вектор правой части преобразуется следующим образом: Ьс =ТЬпс. При этом матрица Ас и вектор

правой части Ьс получаются не прямым умножением на матрицу перехода Т, а с использованием специальных алгоритмов [Ю.Г. Соловейчик]. Для случая использования А - V -постановки на несогласованных сетках в диссертационной работе разработан специальный алгоритм сборки матрицы СЛАУ с одновременным использованием двух матриц перехода: Те - для реберных базисных функций и Т„ - для узловых базисных функций.

Третья глава посвящена подсистеме ЭР-МТЗ программного комплекса GeoEM. Программный комплекс ЭР-МТЗ предназначен для моделирования трехмерных электромагнитных полей в задачах магнитотеллурических зондирований и может быть использован для интерпретации данных и проектирования работ методом магнитотеллурического зондирования. Основным приоритетом при создании программного комплекса ЭР-МТЗ являлась высокая скорость получения решения при соблюдении необходимой точности расчета магнитотеллурических полей. Для достижения этих целей в программном комплексе ЭР-МТЗ реализованы следующие методы и технологии, которые выделяют его из существующих программ конечноэлементного моделирования:

- метод разделения поля на нормальную (поле в горизонтально-слоистой среде) и аномальную (поле от трехмерных включений) составляющие;

- методы автоматического построения оптимизированных трехмерных сеток из параллелепипедов и шестигранников с терминальными узлами;

- метод определения необходимого размера расчетной области;

- построение дискретных аналогов вариационных формулировок для узлового и векторного МКЭ, включая переключение между ними при решении задач на различных частотах для увеличения скорости получения решения;

- эффективные методы решения конечноэлементных СЛАУ с комплексными матрицами, хранящимися в специальных блочных форматах, в том числе реализации этих методов для многоядерных компьютеров;

- удобный интерфейс графического препроцессора, позволяющий пользователю быстро задавать параметры расчетной области;

- пост-обработка конечноэлементного решения для повышения точности выдачи характеристик моделируемых полей и их производных.

Реализация программного комплекса ЭР-МТЗ выполнена на языке программирования С++. Компиляция исполняемого модуля выполнена в среде MS Visual Studio 2005 для операционной системы MS Windows. Структуры данных и методы для работы с ними в основных вычислительных модулях бы-

ли реализованы на базе классов С++. Общая структура программного комплекса представлена на рис. 1. Модули, разработанные автором, выделены жирной рамочкой.

Д2

Д5

Графический — процессор

ДО

Д1

модуль параллельного запуска по гармоникам

Построение ЗО-сетки

Д7

I Д4,.

Д8

Чтение ЗЕ)-сетки и данных препроцессора

Д6

Решение 1Р-задачи

Д9

Построение Т-матрицы для векторного МКЭ

Дп

Д12

Д13

Д14

Генерация портрета матрицы СЛАУ в векторном МКЭ

Д20

Сборка СЛАУ в векторном МКЭ (А-постановка)

Д22

Модуль вычисления локальных матриц л выдачи с конечного элемента в векторном МКЭ

Д15

Построение Т-матрицы для узлового МКЭ

Д16

Д17

Д18

Д19

Генерация портрета матрицы СЛАУ в узловом МКЭ_

Сборка СЛАУ в векторном МКЭ А-У-постановка)

Д21

ИЗ

Сборка СЛАУ в узловом МКЭ

учет краевых условии для векторного МКЭ (А-постановка)

Д25

Модули решения СЛАУ для векторного МКЭ

Д28

Д10

Пост-обработка решения в векторном МКЭ

И4

Д23

Д24

учет краевых условии для векторного МКЭ (А-У-постановка)

Д26

* Модуль вычисления локальных матриц к выдачи с конечного элемента в векторном МКЭ

учет краевых условий для узлового МКЭ

Д27

Модули решения СЛАУ для узлового МКЭ_

Д29

Пост-обработкг решения в узловом МКЭ

И5

программный модуль

• — передача данных

■ использование программных модулей

| | — программные модули, разработанные автором

Рис. 1. Общая структура программного комплекса ЭР-МТЗ

В программном комплексе ЭР-МТЗ реализован алгоритм, позволяющий произвести подбор оптимального размера расчетной области и уменьшить вычислительные затраты при выполнении расчетов трехмерных магнитотеллури-ческих полей.

Идея, на которой базируется данный алгоритм, основывается на явлении так называемого скин-эффекта: изменения глубины проникновения электромагнитного поля с изменением частоты. Таким образом, для расчета поля на низких частотах заведомо требуется задание области большего размера, чем при расчете поля на высоких частотах. Поэтому сначала поле вычисляется на самой низкой частоте, после этого начинается подбор размера "бака", который осуществляется следующим способом. В конечноэлементной сетке ищутся узлы, соседние с узлами, лежащими на границе области. Линейные алгебраические уравнения в конечноэлементной СЛАУ, соответствующие найденным узлам, заменяются на уравнения, обеспечивающие выполнение главного краевого условия. Таким образом, размер "бака" уменьшается и фактическое число неизвестных в задаче сокращается. Далее решается полученная СЛАУ и вычисляются характеристики поля в приемниках (для задач магнитотеллурического зондирования - кажущееся сопротивление р и фаза <р).

Если максимальное отличие в приемниках в области с уменьшенным "баком" и без уменьшения "бака" не превышает заданной относительной погрешности (как правило, это е = 1%), то граница большого "бака" перемещается внутрь области и происходит следующая попытка ее перемещения, если нет, то происходит откат к предыдущей границе "бака" и производится расчет поля на следующей (более высокой частоте).

Проведена верификация решений трехмерных задач магнитотеллуриче-ских зондирований, полученных с использованием программного комплекса, на моделях меньшей размерности и сравнением с результатами международного проекта СОМЕММ1.

В четвертой главе приведены результаты вычислительных экспериментов по моделированию трехмерных геоэлектромагнитных полей. Приводятся результаты исследования временных затрат, требуемых для решения конечно-элементных СЛАУ, получающихся при использовании узлового МКЭ и векторного МКЭ в различных постановках, в зависимости от частоты поля и электрической проводимости среды. Проводится сравнение вычислительных затрат, требуемых для решения конечноэлементных СЛАУ, при использовании различных методов решения СЛАУ, различных предобусловливателей, различных форматов хранения матриц СЛАУ и при запуске решения СЛАУ в многоядерных вычислительных системах. Приведены результаты применения разработанных методов на практических задачах трехмерных магнитотеллурических зондирований на площадях Восточной Сибири и Дальнего Востока.

Приведем результаты вычислительных экспериментов на геоэлектрической модели, представляющей собой 5-слойную среду с удельным сопротивлением слоев в диапазоне от 18 до 20000 Ом-м и 9-ю трехмерными включениями с удельным сопротивлением этих ЗЭ-включений в диапазоне от 2 до 2000 Ом-м.

Построенная для аппроксимации рассматриваемой трехмерной геоэлектрической модели конечноэлементная параллелепипеидальная сетка с терминальными узлами содержала 21633 узла, 18533 элемента и 65569 ребер.

В таблице 1 представлено время решения конечноэлементных СЛАУ, полученных в результате конечноэлементной аппроксимации рассматриваемой задачи, при их последовательном решении на одном ядре компьютера Intel(R) Соге(ТМ) 2 Quad CPU Q6600 2.40 GHz. В качестве критерия останова решателей было взято уменьшение относительной невязки в 104 раз. Для предобуслов-ливания СЛАУ использовались главная диагональ матрицы, состоящая из комплексных чисел.

Время решения конечноэлементных СЛАУ

Таблица 1

У, Гц Время

COCR, векторный МКЭ GMRES(50), узловой МКЭ МНБР, узловой МКЭ

3.81Е-04 19 мин 58 с 2 мин 5 с Зс

1.53Е-03 17 мин 38 с 2 мин 3 с 2с

6.10Е-03 11 мин 27 с 2 мин 7 с 3 с

2.44Е-02 7 мин 0 с 2 мин 8 с 5с

9.77Е-02 4 мин 0 с 1 мин 46 с 8с

0.390625 2 мин 14 с 3 мин 56 с 40 с

1.5625 53 с 11 мин 21 с не сходится

6.25 25 с 9 мин 6 с не сходится

25 13 с 8 мин 58 с не сходится

100 8с 6 мин 58 с не сходится

500 5с 4 мин 36 с не сходится

Таким образом, при использовании векторного МКЭ число итераций, требуемых для решения СЛАУ, значительно возрастало при уменьшении частоты поля, при этом наименьшее время решения СЛАУ для диапазона высоких частот наблюдалось при использовании метода СОСЯ. Наименьшее время решения СЛАУ на низких частотах наблюдалось при использовании постановки для узлового МКЭ и метода МНБР в качестве метода решения полученной СЛАУ.

Проведем результаты проведенных экспериментов по сравнению А -постановки, А — V -постановки и метода деревьев-кодеревьев в векторном МКЭ. Геоэлектрическая модель представляла собой 4-слойную среду с 9 трехмерными объектами с проводимостью от 1 од 1400 Ом-м (эта геоэлектрическая модель была получена на промежуточном этапе интерпретации данных МТЗ, снятых вдоль одного из региональных профилей). Для данной задачи была проведена серия расчетов аномального поля при использовании следующих ВМКЭ-постановок:

1. Использовалась вариационная постановка (2) для вектор-потенциала Аа без использования технологии деревьев-кодеревьев и без введения скалярного потенциала ( А -постановка).

2. Использовалась вариационная постановка (3)-(4) для вектор-потенциала Аа с

введением скалярного потенциала Vа (А -К-постановка).

3. Для А -постановки использовалась технология деревьев-кодеревьев, причем в дерево включались только ребра, находящиеся в подобласти с а = 0. Неизвестные, соответствующие ребрам дерева, исключались.

4. Для А -постановки использовалась технология деревьев-кодеревьев, причем в дерево были включены ребра, находящиеся как в подобласти с о - 0, так и в подобластях с а * 0. Неизвестные, связанные с ребрами дерева исключались, а в узлах сетки, лежащих в подобласти с а * 0, был введен скалярный

потенциал V".

Время решения СЛАУ, соответствующих различным частотам, показано на графике на рисунке 2.

Рис 2. Время решения конечноэлементных СЛАУ в зависимости от частоты поля при использовании ^-постановки, Л-К-постановки и при построении У-дерева во всей расчетной области и 2-дерева в подобласти с а = О

В верхней части частотного диапазона 0.3-500 Гц использование А-постановки и подхода, основанного на построении остовного дерева в подобласти с а = 0, дает сравнимый результат. При этом подходы, основанные на использовании Л-К-постановки и построении остовного дерева во всей расчетной области, требуют большего времени решения СЛАУ. В нижней части частотного диапазона 0.00034-0.03 Гц наименьшее время решения конечноэлементных СЛАУ было получено при использовании Л-К-постановки. При этом для частот

ниже 0.005 Гц время решения СЛАУ для /4-К-постановки с уменьшением частоты не увеличивалось. На самых низких частотах 0.00034-0.002 Гц применение подхода с построением остовного дерева во всей расчетной области позволило решить конечноэлементные СЛАУ за меньшее время, чем при использовании Л-постановки, однако использование A-V-постановки дало еще лучший результат. Так, время решения СЛАУ методом COCR с диагональным предобуслов-ливанием для частоты 0.00034 Гц при построении У-дерева во всей расчетной области составило 19 мин 37 с (5448 итераций), при использовании А-постановки - 2 ч 36 мин (44540 итераций), при использовании Л-К-постановки - 10 мин 54 с (1880 итераций).

Для распараллеливания задач частотных зондирований был реализован подход, основанный на параллельном запуске подзадач, соответствующим разным частотам, на отдельных узлах вычислительной системы. Автором выполнена параллельная реализация программы для систем с распределенной памятью на основе технологии Message Passing Interface (MPI). Для систем с общей памятью также выполнено распараллеливание методов решения СЛАУ на ОрепМР.

Приведем результаты вычислительных экспериментов на геоэлектрической модели, полученной на одном из этапов интерпретации данных МТЗ в районе параметрической скважины в районе Верхне-Илимпейской площади Восточной Сибири. В геоэлектрическую модель включены массивные глубинные аномалии и множество относительно небольших аномалий, расположенных на небольших глубинах (всего 36 объектов). Вид геоэлектрической модели показан на рисунке 3. Конечноэлементная сетка содержит 588007 узлов и 1781314 ребер. Время решения задач для 18 частот в диапазоне от 0.000584 до 316.2 Гц и двух направлений тока (всего 36 задач) при использовании различных методов приведено в табл.2 (на 1 ядре и 10 ядрах). Распараллеливание выполнялось по гармоникам.

Рис 3. 3D-bh4 геоэлектрической модели с 36 объектами

Таблица 2

Время решения задачи на всех частотах

Метод решения Время (1 ядро) Время (10 ядер)

Векторный МКЭ, СОСЯ, А-постановка 439 ч 16 мин 61 ч 44 мин

Векторный МКЭ, СОСЯ, А-У-постановка 46 ч 20 мин 4 ч 45 мин

Комбинирование методов МНБР и СОСЯ (А-постановка) 105 ч 16 мин 17 ч 34 мин

Комбинирование методов МНБР и СОСЯ (А-У-постановка) 28 ч 26 мин 3 ч 37 мин

При использовании метода ОМЯЕБ решение полученных СЛАУ удалось получить только на 4 нижних частотах в диапазоне от 0.000584 до 0.004642 Гц. При этом суммарное время счета на 1 ядре составило более 63 часов. Для более высоких частот метод ОМГШБ не сходился. Методы ВСО и ЛОС не сходились для всех 18 гармоник.

Таким образом, наиболее выгодным по вычислительным затратам оказалось комбинирование узлового МКЭ с решателем МНБР на низких частотах и векторного МКЭ в А-У-постановке с решателем СОСЯ на средних и высоких частотах. При использовании А-У-постановки вместо А-постановки наблюдалось более высокое ускорение за счет распараллеливания задач, соответствующих разным гармоникам, поскольку время решения задач на низких и средних частотах при использовании А-У-постановки на порядок меньше, чем при использовании А постановки.

Эффективность разработанного метода оптимизации размера "бака" была исследована на двух задачах моделирования трехмерных магнитотеллуриче-ских полей. В первой задаче геоэлектрическая модель состояла из 4-слойной среды с сопротивлением слоев от 100 до 10000 Ом-м и 9 трехмерных объектов с сопротивлением от 5 до 10000 Ом-м. Вторая геоэлектрическая модель состояла из 4-слойной среды с сопротивлением слоев от 30 до 800 Ом-м и 13 объектов с сопротивлением от 2 до 70 Ом-м. Начальный размер расчетной области в

обоих случаях составлял [-500...500]3 км3. Суммарное время решения задач на всех частотах при токе, направленном вдоль оси X, составило для первой модели 1 ч 55 мин без применения алгоритма оптимизации размера расчетной области и 31.5 мин с оптимизацией размера расчетной области. Для второй модели — 4 ч 3 мин и 2 ч 37 мин соответственно. Таким образом, можно заключить, что эффективность предложенного алгоритма подбора "бака" оказывается больше при расчетах в средах с высоким удельным сопротивлением. Найденный размер расчетной области составил при этом

[160608,187652] х [-188703,217734] х[-169844,127761] м3 для первой модели и

[-47975,46149] х [-44413,54004] х [-51872,15697] м 3 для второй модели.

Далее сравним методы решения конечноэлементных СЛАУ в трехмерных задачах аэроэлектроразведки на примере геоэлектрической модели, состоящей

из слоя с толщиной 500 м и сопротивлением 200 Ом-м, фундамента с сопротивлением 1000 Ом-м и четырех трехмерных включений с сопротивлением от 1 до 10 Ом-м, лежащих на глубине от 30 до 500 м. Источником служила электрическая петля радиусом 15 м на высоте 50 м. Приемник находился на той же высоте на расстоянии 50 м позади генераторной конструкции. Конечноэлементная сетка содержала 14858 несогласованных конечных элементов (в виде параллелепипедов). Вычислительные затраты, требуемые для решения задач, соответствующих различным гармоникам, при использовании узлового МКЭ и векторного МКЭ с использованием различных методов решения СЛАУ приведены в таблице 3.

Таблица 3

Вычислительные затраты требуемые для решения СЛАУ методами

ОМКЕ5(Ю) (узловой МКЭ) и СОСО, СОСЯ (векторный МКЭ)

частота, Гц ОМКЕБ(10) СОСО СОСИ.

итерации время итерации время итерации время

130 2585 9 мин 44 с 3463 1 мин 2585 44 с

2080 52 11 с 915 15с 809 14 с

8320 14 3 с 522 9с 436 7с

На основании этих данных можно заключить, что метод ОМКЕБ для данной задачи по прежнему проигрывает по времени решения конечноэлементных СЛАУ методам СОСО и СОСЯ, а метод СОСЯ позволяет получить решение конечноэлементной СЛАУ за несколько меньшее время, чем метод СОСО.

Заключение

Основные результаты работы заключаются в следующем.

1. Разработаны методы моделирования трехмерных гармонических по времени электромагнитных полей с использованием узлового и векторного МКЭ с выделением части поля на сетках с терминальными узлами.

2. Показано, что использование метода деревьев-кодеревьев при моделировании трехмерных магнитотеллурических полей векторным МКЭ позволяет снизить вычислительные затраты при расчетах в нижней части диапазона частот при условии включения в кодерево ребер конечноэлементной сетки, совпадающих с направлением стороннего тока. Однако гораздо большее снижение вычислительных затрат при использовании векторного МКЭ в нижнем и в среднем диапазоне частот было достигнуто с применением математической постановки с дополнительным введением скалярного потенциала. В верхнем диапазоне частот наименьшие затраты при решении конечноэлементных СЛАУ наблюдались для математической постановки с использованием одного вектор-потенциала.

3. Разработан алгоритм решения конечноэлементных СЛАУ при расчетах трехмерных электромагнитных полей на низких частотах — метод минимизации невязки с использованием блочной релаксации (МНБР). Показаны его преимущества при расчетах трехмерных магнитотеллурических полей на низких час-

тотах по сравнению с другими методами решения СЛАУ и в сравнении с использованием математической постановки для векторного МКЭ.

4. Показано, что при моделировании трехмерных электромагнитных полей в задачах магнитотеллурических зондирований комбинирование узлового МКЭ с решателем МНБР и векторного МКЭ в постановке с совместным использованием векторного и скалярного потенциалов с решателем COCR позволяют добиться наименьшего времени решения задачи.

5. Разработана и реализована вычислительная часть подсистемы ЭР-МТЗ программного комплекса GeoEM, предназначенная для моделирования трехмерных магнитотеллурических полей и включающая в себя модули генерации и решения конечноэлементных СЛАУ для различных постановок с использованием скалярных и векторных базисных функций.

6. Программный комплекс ЭР-МТЗ применялся для интерпретации данных на площадях Восточной Сибири при нефтепоисковых исследованиях и для интерпретации данных, полученных при геофизических исследованиях на региональных профилях Дальнего Востока.

Основные результаты диссертации опубликованы в следующих работах:

1. Домников П.А., C.B. Киреева, М.Г. Персова, Ю.Г. Соловейчик. Конечно-элементное моделирование трехмерных магнитотеллурических полей с применением технологии деревьев-кодеревьев и постановки с совместным использованием векторного и скалярного потенциалов // Научный вестник НГТУ. - 2011. - №3. - С. 43-52.

2. Домников П.А., М.Г. Персова, Ю.Г. Соловейчик, Д.В. Вагин. Моделирование трехмерных магнитотеллурических полей векторным методом конечных элементов и возможности распараллеливания на процессорах с общей памятью И Научный вестник НГТУ. - 2010. - № 3 (40). - С. 87-96.

3. Домников П.А. Распараллеливание задач конечноэлементного моделирования трехмерных магнитотеллурических полей // Информационные технологии. - 2010. - №3. - С. 69-72.

4. Домников П.А., М.Г. Персова, Ю.Г. Соловейчик, Ю.В. Тракимус. Моделирование трехмерных нестационарных электромагнитных полей вертикальной электрической линии и оценка возможности ее использования при поиске залежей углеводородов в условиях шельфовой зоны // Научный вестник НГТУ.-2010.-№3(40).-С. 75-86.

5. Персова М.Г., Ю.Г. Соловейчик, Д.В. Вагин, П.А. Домников. Сравнение различных подходов к численному моделированию трехмерных полей вызванной поляризации // Доклады АН ВШ. - 2011. — №2. — С. 126-139.

6. Персова М.Г., Ю.Г. Соловейчик, П.А. Домников, А.Г. Задорожный. Сравнение возможностей векторного и узлового МКЭ при решении трехмерных задач магнитотеллурических зондирований И Научный вестник НГТУ. -2011.-№2(43).-С. 83-96.

7. Персова М.Г., Ю.Г. Соловейчик, М.В. Абрамов, Д.В. Вагин, П.А. Домников. Программные комплексы конечноэлементного моделирования электромагнитных полей в технических устройствах и в задачах геоэлектрики //

Вестник НГУ. Серия: Информационные технологии. - 2012. - Том 10. -№2.-С. 34-43.

8. Домников П.А. Метод решения систем уравнений, возникающих при ко-нечноэлементной аппроксимации гармонических по времени электромагнитных полей // Сборник научных трудов Hl ТУ. - 2009. №2(56). - С. 41-46.

9. Домников П.А. Моделирование трехмерных электромагнитных полей для магнитотеллурических зондирований узловым и векторным МКЭ // Сборник научных трудов НГТУ.-2008.-№ 1(51).-С. 55-60.

10.Персова М.Г., Ю.Г. Соловейчик, Д.В. Вагин, П.А. Домников, М.Г. Токарева, Г.М. Тригубович. О подходе к решению прямых и обратных трехмерных задач вызванной поляризации [Электронный ресурс] // Материалы международной научно-практической конференции Геобайкал 2012. - г. Иркутск.

- 20-24 августа 2012. - EAGE. - 5 С. - Режим доступа: http:// earthdoc.eage.org/detail.php?pubid=61526.

11 .М.Г. Персова, Ю.Г. Соловейчик, Д.В. Вагин, П.А. Домников, М.Г. Токарева. Реализация решения научных и производственных задач электроразведки в распределенных вычислительных системах // Материалы международной научно-практической конференции Геобайкал 2012. - г. Иркутск. - 2024 августа 2012. - EAGE. - 5 С. - Режим доступа: http://earthdoc.eage.org/ detail.php?pubid=61527.

12.М.Г. Персова, Ю.Г. Соловейчик, Д.В. Вагин, П.А. Домников, М.Г. Токарева Особенности 30-моделирования при решении задач морской электроразведки // Материалы международной научно-практической конференции Геобайкал 2012. - г. Иркутск. - 20-24 августа 2012. - EAGE. - 5 С. - Режим доступа: http://earthdoc.eage.org/detail.php?pubid=61528.

13.Домников П.А. О решении задач конечноэлементного 30-моделирования низкочастотных магнитотеллурических полей в изотропной и анизотропной среде // Математическое и компьютерное моделирование естественнонаучных и социальных проблем: сборник статей IV международной научно-технической конференции молодых специалистов, аспирантов и студентов.

- Пенза: Приволжский Дом знаний, 2010. - С. 184-185.

14.Домннков П.А. Моделирование трехмерных магнитотеллурических полей векторным методом конечных элементов // Актуальные проблемы электронного приборостроения. АПЭП-2008.: материалы IX международной конференции, Новосибирск, 24-26 сентября, 2008. — Новосибирск: Изд-во НГТУ, 2008. - Том 6. - С. 93-98.

15.Домников П.А., C.B. Киреева, И.Е. Аврунева. О методах решения конечно-элементных СЛАУ при ЗО-моделировании низкочастотных магнитотеллурических полей // Актуальные проблемы электронного приборостроения. АПЭП-2010.: материалы X международной конференции, Новосибирск, 2224 сентября, 2010. - Новосибирск: Изд-во НГТУ, 2010. - Том 6. - С. 118-123.

16.Домников П.А., C.B. Киреева. О решении задач электромагнетизма векторным методом конечных элементов с вырожденной матрицей жесткости // Российская научно-техническая конференция "Информатика и проблемы

телекоммуникаций": материалы конференции, Новосибирск, 27-28 апреля 2010. - Новосибирск: Изд-во СибГУТИ, 2010. - Том 1. - С. 55-57.

17.ДОМНИКОВ П.А. Об одном методе решения систем конечноэлементных уравнений, возникающих в задачах моделирования низкочастотных электромагнитных полей // Наука. Технологии. Инновации.: материалы всероссийской научной конференции молодых ученых, Новосибирск, 4-7 декабря, 2008. - Новосибирск: Изд-во НГТУ, 2008. - Часть 1. - С. 72-74.

18.Домников П.А. Об особенностях конечноэлементных постановок для трехмерных гармонических по времени электромагнитных полей в анизотропных средах // Наука. Технологии. Инновации.: материалы всероссийской научной конференции молодых ученых, Новосибирск, 4-5 декабря, 2009. -Новосибирск: Изд-во НГТУ, 2009. - Часть 1. - С. 87-89.

19.Домников П.А. Организация параллельных вычислений в задачах конечно-элементного моделирования трехмерных магнитотеллурических полей // Технологии Microsoft в теории и практике программирования: сборник трудов VI Всероссийской научно-практической конференции студентов, аспирантов и молодых ученых. Томск, 17-18 марта 2009. - Томск: Изд-во Томского политехнического университета, 2009. - С. 50-53.

20.Домн11ков П.А. Решение трехмерных задач магнитотеллурического зондирования векторным методом конечных элементов // Российская научно-техническая конференция "Информатика и проблемы телекоммуникаций": материалы конференции, Новосибирск, 26-28 апреля 2007. - Новосибирск: Изд-во СибГУТИ, 2007. -Том 1. - С. 119-120.

21. Домников П.А. Сравнение векторного и узлового метода конечных элементов для моделирования трехмерных электромагнитных полей в задачах магнитотеллурического зондирования // Российская научно-техническая конференция "Информатика и проблемы телекоммуникаций": материалы конференции, Новосибирск, 24-25 апреля 2008. - Новосибирск: Изд-во СибГУТИ, 2008. - Том 1. - С. 45-48.

22.Домников П.А.Распараллеливание процесса решения задач конечноэле-ментного моделирования трехмерных гармонических по времени геоэлектромагнитных полей // Наука. Технологии. Инновации.: материалы всероссийской научной конференции молодых ученых, Новосибирск, 3-5 декабря, 2010.-Часть 1.-С. 61-63.

23.ПерсоваМ.Г., Ю.Г. Соловейчик, Г.М. Тригубович, Д.В. Вагин, П.А. Домников Программное обеспечение для моделирования трехмерных электромагнитных полей в различных технологиях электроразведки // Материалы Всероссийской школы-семинара им. М.Н. Бердичевского и Л.Л. Ваньяна по электромагнитным зондированиям Земли, 16-21 мая 2011. -ЭМЗ-2011. В двух книгах. Книга 2. - Спб.: СПбГУ, 2011. - С. 467-471.

24.Киреева C.B., П.А. Домников. Решение трехмерной задачи магнитотеллурического зондирования при совместном использовании векторных и скалярных базисных функций // Российская научно-техническая конференция "Информатика и проблемы телекоммуникаций": материалы конференции, Новосибирск, 21-22 апреля 2011. - СибГУТИ, 2011. - Том 1. - С. 81-83.

25.Персова М.Г., Ю.Г. Соловейчик, П.А. Домников, C.B. Киреева. 3D-моделирование магнитотеллурических полей с использованием распараллеливания // Актуальные проблемы электронного приборостроения. АПЭП-2010.: материалы X международной конференции, Новосибирск, 22-24 сентября, 2010. - Новосибирск: Изд-во НГТУ, 2010. - Том 6. - С. 144-149.

26.Персова М.Г., Ю.Г. Соловейчик, П.А. Домников, Е.И. Симон. Компьютерное моделирование трехмерных геоэлектромагнитных полей, возбуждаемых контролируемыми источниками, в средах с анизотропной проводимостью // Актуальные проблемы электронного приборостроения. АПЭП-2010.: материалы X международной конференции, Новосибирск, 22-24 сентября, 2010. - Новосибирск: Изд-во НГТУ, 2010. - Том 6. - С. 134-138.

27.Персова М.Г., Ю.Г. Соловейчик, П.А. Домников, Е.И. Симон. О томографическом подходе к интерпретации данных электромагнитных зондирований в трехмерных средах // Актуальные проблемы электронного приборостроения. АПЭП-2010.: материалы X международной конференции, Новосибирск, 22-24 сентября, 2010. - Новосибирск: Изд-во НГТУ, 2010. - Том 6. -С. 150-154.

28.Соловейчик Ю.Г., Персова М.Г., Абрамов М.В., Вагин Д.В., Домннков П.А. HEL-MAR - ЗО-моделирование геоэлектромагнитных полей горизонтальной электрической линии для решения задач морской электроразведки // Свидетельство о государственной регистрации программы для ЭВМ № 2012618323. М.: Федеральная служба по интеллектуальной собственности (Роспатент). - 2012.

29.Соловейчик Ю.Г., Персова М.Г., Абрамов М.В., Вагин Д.В., Домннков П.А., Токарева М.Г. GeoEM (GeoElectroMagnetics) // Свидетельство о государственной регистрации программы для ЭВМ № 2012614504. М.: Федеральная служба по интеллектуальной собственности (Роспатент). - 2012.

Отпечатано в типографии Новосибирского государственного технического университета 630092, г. Новосибирск, пр. К. Маркса, 20,

тел./факс (383) 346-08-57 формат 60 X 84/16 объем 1.5 п.л. тираж 100 экз. Заказ № 1544 подписано в печать 9.11.2012 г

Введение

Глава 1. Математические модели гармонических электромагнитных полей с использованием МКЭ.

1.1. Математическое моделирование электромагнитных полей в задачах геоэлектромагнетизма.

1.2. Задачи магнитотеллурических зондирований.

1.3. Математические модели с выделением поля вмещающей среды.

1.4. Эквивалентные вариационные постановки.

1.5. Выводы.

Глава 2. Построение дискретных аналогов и разработка методов решения конечноэлементных СЛАУ.

2.1. Построение дискретных аналогов в узловом МКЭ.

2.2. Методы решения конечноэлементных СЛАУ, возникающих при использовании узлового МКЭ.

2.3. Метод минимизации невязки с использованием блочной релаксации.

2.4. Построение дискретного аналога для постановки с одним вектор-потенциалом в векторном МКЭ.

2.5. Построение дискретного аналога для вариационной постановки с введением градиента скалярного потенциала.

2.6. Построение дискретных аналогов с применением метода деревьевкод еревьев

2.7. Методы решения конечноэлементных СЛАУ, возникающих при использовании векторного МКЭ.

2.8. Выводы.

Глава 3. Описание программного комплекса ЭР-МТЗ.

3.1. Назначение и общая структура программного комплекса.

3.2. Описание основных программных модулей.

3.3. Описание модулей решения СЛАУ.

3.4. Модуль оптимизации размеров расчетной области.

3.5. Верификация программного комплекса.

3.6. Выводы.

Глава 4. Результаты вычислительных экспериментов.

4.1. Исследование времени решения СЛАУ в зависимости от частоты поля и проводимости среды.

4.2. Сравнение производительности разработанных подпрограмм с библиотекой математических функций Intel MKL.

4.3. Компьютерная реализация методов решения СЛАУ в многоядерных системах.

4.4. Применение алгоритма подбора размера расчетной области.

4.5. Применение метода деревьев-кодеревьев и постановки с совместным использованием скалярного и векторного потенциалов в векторном МКЭ.

4.6. Применение разработанных методов для решения практических задач магнитотеллурических зондирований на площадях Восточной Сибири и

Дальнего Востока.

4.7. Решение конечноэлементных СЛАУ в трехмерных задачах с источником в виде горизонтальной электрической линии.

4.8. Решение конечноэлементных СЛАУ в трехмерных задачах аэроэлектроразведки.

4.9. Выводы.

Задачи моделирования трехмерных геоэлектромагнитных полей в неоднородных средах относятся к одним из наиболее вычислительно трудоемких. Уменьшение вычислительных затрат и повышение точности решения этих задач, несомненно, является актуальной задачей, поскольку от этого напрямую зависит качество, а часто и сама возможность адекватного восстановления геоэлектрического строения среды.

Основными методами, применяемыми для численного моделирования трехмерных геоэлектромагнитных полей, являются метод интегральных уравнений (МИУ), метод конечных разностей (МКР) и метод конечных элементов (МКЭ). При этом известно, что МИУ, в разработку которого существенный вклад внесли [5, 6, 17, 18, 49, 70, 76, 77, 108, 109, 131 и др.], является вычислительно эффективным при расчетах полей в средах с небольшим количеством трехмерных объектов, с ростом же количества трехмерных неоднородностей (что более соответствует реальным условиям) этот метод начинает существенно уступать по вычислительной эффективности таким сеточным методам, как МКР (для задач геоэлектрики он разрабатывался [19, 39, 48, 68, 72, 107 и др.]) и МКЭ, в развитие которого для решения задач геоэлектромагнетизма значимый вклад внесли [14, 23, 42, 45, 64, 118, 128, 130 и др.]. МКР и МКЭ основаны на аппроксимациях дифференциальных уравнений в частных производных. При этом по вычислительной эффективности МКР, в свою очередь, довольно сильно уступает МКЭ из-за необходимости использования в МКР прямоугольных регулярных сеток - для получения решения требуемой точности это приводит к дискретным аналогам гораздо большей размерности.

В принципе, МКР и МКЭ уже достаточно давно применяются для решения трехмерных задач электроразведки. Основным преимуществом этих методов перед МИУ является разреженность матриц систем уравнений, получаемых в результате аппроксимаций соответствующих трехмерных задач. Но вместе с тем МКР и МКЭ требуют включения в расчетную область помимо самих трехмерных объектов довольно большого пространства вокруг них, и при этом для достижения хорошей точности необходимы достаточно подробные сетки, причем мелкие ячейки нужны не только в ЗБ-объектах, но и в окружающем их пространстве, особенно вблизи источников поля и приемников. В результате, при решении задач с контролируемыми источниками, если ЗБ-объекты в них дают относительно слабые отклики, МКР и МКЭ требуют очень высоких вычислительных затрат на получение решения с приемлемой точностью.

Однако разрабатываемые в последние годы модификации МКЭ с включением в него возможности выделения поля простой структуры (что является стандартным для МИГУ и не использовалось до самого последнего времени в МКЭ) - так называемого нормального поля или поля вмещающей среды - существенно расширили его возможности при решении трехмерных задач геоэлектрики. Присущие МКЭ возможности учета геометрии расчетной области любого уровня сложности и его гибкость при построении сеток в сочетании с технологиями, использующими выделение нормального поля, позволяют сделать этот метод мощным и удобным инструментом при решении трехмерных задач геоэлектромагнетизма.

При использовании МКЭ для моделирования трехмерных геоэлектромагнитных полей основная вычислительная нагрузка, как правило, приходится на решение систем линейных алгебраических уравнений (СЛАУ), получаемых в результате конечноэлементной аппроксимации исходной краевой задачи.

При решении трехмерных задач, особенно векторных, конечноэлемент-ные аппроксимации приводят к СЛАУ с разреженными матрицами очень большого размера с довольно большой шириной ленты, поэтому для решения таких СЛАУ применяют, в основном, итерационные методы, которые рассмотрены, например, в работах [7, 13, 55, 67, 88, 112, 115 и др.]. При этом скорость сходимости различных итерационных методов может сильно (в десятки, а иногда и в сотни раз) различаться, причем при решении конкретной практической задачи, как правило, заранее неизвестно, какой из методов решения может оказаться наиболее быстрым. Поэтому исследование различных методов решения СЛАУ, получаемых в результате конечноэлементной аппроксимации трехмерных задач, также является актуальной задачей.

При моделировании трехмерных векторных полей в методе конечных элементов наиболее широко распространенными являются следующие два подхода. Первый из них (узловой или скалярный МКЭ) состоит в том, что каждая из компонент искомого поля ищется в виде линейной комбинации скалярных базисных функций (которые обычно связывают с узлами конечноэлементной сетки). Второй подход (векторный МКЭ) позволяет строить математические модели непосредственно в терминах векторных переменных (используются векторные базисные функции). В условиях, когда токами смещения можно пренебречь (моделирование электромагнитных полей на не очень высоких частотах), можно использовать постановки, как для векторных, так и для скалярных базисных функций, и в зависимости от решаемой задачи преимущество может иметь как одна, так и другая постановка.

Отметим, что даже при использовании технологии выделения поля вычислительные затраты, требуемые для расчета трехмерных полей в сложных средах, содержащих много трехмерных включений, могут оставаться довольно большими. Дополнительный эффект сокращения вычислительных затрат может быть получен за счет оптимизации сеток с обоснованным изменением шагов и устранением так называемых "лишних" узлов из регулярных сеток. В этом случае в сетках могут появляться терминальные узлы, и необходима разработка специальных подходов к построению дискретных аналогов. Кроме того, вычислительные затраты при решении задач геоэлектромагнетизма с использованием МКЭ довольно сильно зависят от размера расчетной области. Ее границы должны быть удалены от мест регистрации сигнала настолько, чтобы краевые условия практически не влияли на получаемые в приемниках значения, поэтому важной является задача определения удаленных границ расчетной области (границ так называемого большого объема или "бака"). Реализация возможности определения минимально допустимого размера "бака" позволяет дополнительно сократить размер соответствующих дискретных аналогов и, соответственно, уменьшить вычислительные затраты.

В целом, все рассмотренные проблемы, связанные со снижением вычислительных затрат при конечноэлементном моделировании трехмерных геоэлектромагнитных полей, решению которых посвящена данная диссертационная работа, являются отдельными аспектами основной научной проблемы создания эффективных методов численного моделирования трехмерных гармонических по времени электромагнитных полей.

Цели и задачи исследования

1. Построение эффективных конечноэлементных аппроксимаций для решения задачи моделирования трехмерных гармонических по времени электромагнитных полей. Для достижения данной цели сформулированы и решены следующие задачи:

- разработка и программная реализация конечноэлементных вычислительных схем моделирования гармонических по времени электромагнитных полей на трехмерных несогласованных сетках на основе узлового и векторного МКЭ с использованием выделения основной части поля;

- построение дискретных аналогов при использовании метода деревьев-кодеревьев и совместного использования скалярного и векторного потенциалов в векторном МКЭ на трехмерных сетках; реализация разработанных вычислительных схем в программном комплексе.

2. Разработка и реализация эффективных процедур решения конечноэлементных СЛАУ, возникающих при конечноэлементной аппроксимации трехмерных гармонических по времени электромагнитных полей. Для достижения данной цели сформулированы и решены следующие задачи:

- экспериментальное исследование существующих методов решения СЛАУ, возникающих при использовании различных конечноэлементных постановок для моделирования гармонических электромагнитных полей;

- разработка методов решения СЛАУ, получаемых при использовании узлового МКЭ для моделирования трехмерных гармонических электромагнитных полей;

- программная реализация методов решения конечноэлементных СЛАУ, получаемых в результате аппроксимации трехмерных гармонических по времени электромагнитных полей; сравнение производительности разработанных программ с существующими математическими библиотеками.

На защиту выносятся

1. Математические постановки и схемы конечноэлементного моделирования трехмерных гармонических по времени электромагнитных полей с использованием узлового и векторного МКЭ с выделением части поля на сетках с терминальными узлами.

2. Методы решения СЛАУ, возникающих при моделировании трехмерных электромагнитных полей на высоких частотах с использованием векторного МКЭ и на низких частотах с использованием узлового МКЭ.

3. Программные модули конечноэлементного моделирования трехмерных гармонических геоэлектромагнитных полей и решения СЛАУ в подсистеме ЭР-МТЗ программного комплекса ОеоЕМ.

4. Результаты использования разработанных методов и программного комплекса для решения задач геоэлектромагнетизма.

Достоверность результатов

Вычислительные схемы конечноэлементного моделирования трехмерных электромагнитных полей протестированы путем сравнения с решением задач меньшей размерности, с решениями, полученными другими авторами, а также путем сравнения решений одних и тех же задач с использованием представлений электромагнитного поля через различные потенциалы в векторном и в узловом МКЭ.

Научная новизна

1. Разработана и реализована вычислительная схема для моделирования трехмерных магнитотеллурических полей в математической постановке для векторного МКЭ с использованием векторного и скалярного потенциалов, технологии выделения поля вмещающей горизонтально-слоистой среды и нерегулярных сеток с терминальными узлами.

2. Разработана и реализована вычислительная схема для моделирования трехмерных гармонических геоэлектромагнитных полей с использованием метода деревьев-кодеревьев и технологии выделения части поля.

3. Разработан новый метод решения СЛАУ, возникающих при моделировании трехмерных гармонических по времени электромагнитных полей на низких частотах с использованием узлового МКЭ - метод минимизации невязки с использованием блочной релаксации (МНБР). Практическая значимость

Предлагаемые в данной работе способы построения дискретных аналогов и методы решения СЛАУ реализованы в подсистемах ЭР-МТЗ и ЭРИАН программного комплекса ОеоЕМ, который использовался при проектировании поисково-разведочных работ и интерпретации практических данных. Публикации

По результатам выполненных исследований опубликовано 29 работ, в том числе 7 статей в журналах, входящих в перечень изданий, рекомендуемых Высшей аттестационной комиссией для защиты диссертаций, 2 статьи в сборнике научных трудов, 18 публикаций в материалах конференций, 2 свидетельства о государственной регистрации программы для ЭВМ. Апробация работы

Основные результаты работы были представлены и докладывались на: международной научно-практической конференции "Геобайкал 2012" (Иркутск, 2012); всероссийской школе-семинаре им. М.Н. Бердичевского и Л.Л. Ваньяна по электромагнитным зондированиям Земли, ЭМЗ-2011 (Санкт-Петербург, 2011); рабочем семинаре "Наукоемкое программное обеспечение" в рамках Ершовской конференции по информатике (РБГП) (Новосибирск, Академгородок, 2011); IX и X международной конференции "Актуальные проблемы электронного приборостроения" АПЭП-2008, АПЭП-2010 (Новосибирск, 2008,

2010); всероссийской конференции по вычислительной математике КВМ-2011 (Новосибирск, Академгородок, 2011); научной конференции Лаборатории НГУ-Интел "Технологии высокопроизводительных вычислений" (Новосибирск, Академгородок, 2009, 2010); IV международной научно-технической конференции молодых специалистов, аспирантов и студентов "Математическое и компьютерное моделирование естественнонаучных и социальных проблем" (Пенза, 2010); всероссийской научной конференции молодых ученых "Наука. Технологии. Инновации" (Новосибирск, 2008, 2009, 2010); российской научно-технической конференции "Информатика и проблемы телекоммуникаций" (Новосибирск, 2007, 2008, 2010, 2011); VI Всероссийской научно-практической конференции студентов, аспирантов и молодых ученых "Технологии Microsoft в теории и практике программирования" (Томск, 2009).

Личный вклад

Разработан метод решения СЛАУ, возникающих при моделировании гармонических по времени трехмерных электромагнитных полей на низких частотах с использованием узлового МКЭ. Разработаны, реализованы и апробированы программные модули конечноэлементного моделирования и решения СЛАУ в программном комплексе GeoEM. Разработан и реализован алгоритм подбора размера расчетной области для решения трехмерных задач магнитотеллурики. Выполнена компьютерная реализация методов решения конечноэлементных СЛАУ в многоядерных системах с общей памятью и реализация запуска задач магнитотеллурических зондирований по частотам в системах с распределенной памятью. Разработаны и исследованы методы построения дискретных аналогов при моделировании трехмерных электромагнитных полей с использованием метода деревьев-кодеревьев.

В совместных публикациях автору принадлежат следующие результаты. В работах [94, 117, 119] автором проведены расчеты трехмерных электромагнитных полей, в работах [92, 98, 120, 122, 132, 133] автором разработаны методы решения конечноэлементных СЛАУ, в работах [116, 123] автором разработаны и реализованы модули распараллеливания расчетов трехмерных электромагнитных полей и выполнены вычислительные эксперименты, в [90, 97, 113] проведены сравнения итерационных методов решения СЛАУ.

Работа выполнялась при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы (автор принимал участие в 9 проектах, в том числе в 6 проектах являлся руководителем) и ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2013 годы» (в 3 проектах). Работа выполнялась в рамках государственных заданий высшим учебным заведениям в части проведения НИР (шифр заявки 8.874.2011), а также была поддержана грантом Фонда содействия отечественной науке в рамках программы «Лучшие аспиранты РАН» в 2010 г.

Структура работы

Диссертационная работа состоит из введения, четырех глав, заключения, списка использованных источников (133 наименования), приложения и содержит 40 рисунков и 52 таблицы. Общий объем работы -191 страница.

ЗАКЛЮЧЕНИЕ