Параллельная реализация математического моделирования процессов электромагнитного каротажа зондовым комплексом ВИКИЗ Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 12, № 6, 2007

ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ПРОЦЕССОВ ЭЛЕКТРОМАГНИТНОГО КАРОТАЖА ЗОНДОВЫМ КОМПЛЕКСОМ ВИКИЗ

В. Н. Еремин Научно-производственное предприятие геофизической аппаратуры “Луч”, Новосибирск, Россия e-mail: [email protected]

О. В. Нечаев Институт нефтегазовой геологии и геофизики им. А.А.Трофимука СО РАН, Новосибирск, Россия e-mail: [email protected]

Ш. ХАБЕРХАУЭР

NEC — High Performance Computing Europe GmbH, Stuttgart, Germany

e-mail: [email protected]

Н. Ю. ШокинА High Performance Computing Center Stuttgart, Stuttgart, Germany Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

Э. П. Шурина

Новосибирский государственный технический университет, Россия

e-mail: [email protected]

The results of the joint project of the Institute of Computational Technologies of Siberian Branch of Russian Academy of Sciences (ICT SB RAS, Novosibirsk, Russia) and the High Performance Computing Center Stuttgart (HLRS, Stuttgart, Germany) are presented. The project is realized within the framework of the activities of the German-Russian Center for Computational Technologies and High Performance Computing (http://www.grc-hpc.de). The parallel realization of the problem on mathematical modelling of three-dimensional electromagnetic fields in inhomogeneous domains in frequency domain is presented. The computational schemes are investigated in domains,

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

which are inhomogeneous with respect to physical (electroconductivity) and geometrical (different scales of individual fragments of a modelling domain) characteristics, in the presence of field sources changing geometrical position with respect to interfragmentary boundaries, which separate individual subdomains. The three-dimensional vector Helmholtz equation is approximated on inhomogeneous tetrahedral grid with local grid refinements using the vector finite element method on the vector basis functions of high order. The basis functions (edge-elements) for calculating the electric field guarantee continuity of tangential components of the field on interelement and interfragmentary boundaries. For exact calculation of a jump of electric field normal component the special procedure for taking into account a weak divergence on basis functions of high order is realized. The problem on high-frequency induction logging isoparametric soundings (the VIKIZ-method in Russia) in a domain with complicated geometry is solved. The parallelization is done by formulating of n independent problems, where n is a number of probe positions, and associating each of these problems with a separate processor of computing system. The calculation results are presented, obtained on NEC SX-8 cluster.

Введение

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

Метод ВИКИЗ основан на измерении относительных фазовых характеристик электромагнитных величин, а именно, разности фаз наведенных в измерительных катушках электродвижущих сил. Для реализации глубинности исследования, разрешающей способности и чувствительности к определяемым параметрам аппаратура ВИКИЗ и ее модификации (например, скважинный прибор высокочастотного электромагнитного каротажного зондирования в процессе бурения ВИКПБ [2], разработанный в научнопроизводственном предприятии геофизической аппаратуры “Луч” (http://www.looch. ru/index_en.php)) должны иметь несколько зондов высокочастотного электроманит-ного зондирования разной глубинности.

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

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

Область частот в электромагнитном каротаже находится в диапазоне от несколь-

Рис. 1. Схема зондовой части скважинного прибора

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

Для геофизических приложений характерно наличие материалов, электрическое сопротивление и диэлектрическая проницаемость которых контрастно изменяются. Например, в нефтегазоносных районах удельное сопротивление пластов может колебаться в пределах от 2 до 200 Ом-м, а удельное сопротивление бурового раствора в скважине — от 0.01 до 5 Ом-м. Электрическое сопротивление отдельных фрагментов области моделирования контрастно изменяется при переходе через границы этих подобластей, а при работе с обсадными скважинами разрывным является и коэффициент магнитной проницаемости.

Поэтому математическая модель должна выбираться тщательно, с учетом всех сложностей задачи.

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

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

Рис. 2. Характерная область исследования нефтеносных районов : 1 — скважина, заполненная буровым раствором; 2 — диэлектрический корпус зонда, содержащий генераторную катушку Т и приемные катушки Я\ и Д2; 3 — зона проникновения бурового раствора; 4 — пласт; 5 — пропласток

средах с удельным сопротивлением до 120 Ом-м, что расширяет диапазон определяемых удельных сопротивлений.

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

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

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

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

1. Математическая модель

Поведение гармонического по времени электрического поля описывается векторным уравнением Гельмгольца:

го1 — гс^ Е + к2 Е = —гwJ0, (1)

ц

где ц — магнитная проницаемость; к2 = гша — ш2е; Е — напряженность электрического поля, комплексная векторная функция; ш — циклическая частота; е — диэлектрическая проницаемость; а — электрическая проводимость; г — мнимая единица; Jo — плотность тока источника. Необходимо отметить, что электрические характеристики скважины и околоскважинных областей, вмещающих пород и нефте- или газоносных пропластков таковы, что при работе в частотном диапазоне аппаратуры ВИКИЗ токи проводимости и токи смещения соизмеримы. Это приводит к комплексной величине к2 (волновое число) и комплексной векторной функции напряженности электрического поля Е.

Имеет место закон сохранения заряда:

&у ((а + шє)Е) = 0. Также выполняется следующее условие:

30 = 0.

(2)

(3)

Условия непрерывности должны выполняться на границе Г между подобластями с различными материалами:

[п х Е] = 0,

[п ■ (а + г^є)Е]

0,

(4)

(5)

где п — вектор нормали к поверхности Г.

На границе области дП заданы однородные краевые условия:

п х Е

0.

дП

(6)

г

2. Векторная вариационная постановка

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

— непрерывность тангенциальной компоненты поля Е на межфрагментарных границах между подобластями с различными электрическими свойствами;

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

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

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

уравнений Гельмгольца (1), осложняется большим ядром rot-оператора, содержащим все градиенты скалярных функций, принадлежащих пространству И1(П) [3].

В области П введем гильбертовы пространства [4, 5]

H(grad; П) = Н1(П) = {u G L2(n), grad u G L2(n)3},

H(rot; П) = {u G Ь2(П)3, rot u G L2(n)3},

H(div; П) = {u G Ь2(П)3, div u G L2(П)}

с нормами

|u|rot)n = J u ■ u*dQ + J гot u ■ гot u*dQ,

ПП

u|div П = u ■ u*dQ + div u ■ div u*dQ.

1Шу,П

П П

Скалярное произведение определяется следующим образом:

(и, V) = J и •

П

Пространство Н(го1; П) было введено в работах Л.-С. Nёdёlec, свойства этого пространства приведены в [4, 5]. Базисные функции Nёdёlёc — это векторные, го^конформ-ные базисные функции, которые широко используютя для конечно-элементной аппроксимации векторного уравнения Гельгольца и сохраняют непрерывность тангенциальной компоненты векторного поля (например, поля Е) на межэлементных и межфрагментар-ных границах.

Если трехмерная область П (возможно, неоднородная по физическим свойствам) имеет липшиц-непрерывную границу дП, то для любой функции и £ Н(го1; П) можно определить тангенциальный след и х п на дП как элемент Н-1/2(дП) [6].

Функции из пространства Н(го1; П) с нулевым тангенциальным следом формируют его подпространство:

H0(rot; Q) = |u Є H(rot; Q), u x n Аналогично определяются H0(Q) и Ho(div; Q):

H0(Q) = (u Є H1(Q),u =

дП

H0(div; Q) = ju Є H(div; Q), u ■ n

=0

дП

=0

дП

Пространства Щ(П), Н0(го1;П), Н0^^;П), Ь2(П) и три оператора (V, Ух и V-) формируют комплекс де Рама [7], который в Д3 имеет следующий вид:

Н0 А НоМ; П) А Ho(div; П) А Ь2(П). (7)

В соответствии с (7) для пространства Н0(го1;П) имеет место следующее свойство вложения:

Уф £ Н^П), gгad ф £ Н0(го1;П). (8)

Вариационная постановка [3] для задачи (1), (6) имеет вид:

Векторная вариационная постановка. Для ^ Є Ь2(П)3 найти Е Є Н0(го!;П) такое, что для Vv Є Н0(го!; П):

— го! Е, го! ^ + (к2Е, V) = — і (^о, V). (9)

Вариационная постановка (9) выполняется для всех V Є Н0(го!;П). Если, соглас-

но (8), взять V = gгad ф, ф Є Н0(гоі; П), тогда (9) имеет вид

— го! Е, гotgгadф] + (к2Е, gгadф) = —і (^0, gгadф), Vф Є Н0-1 у

Учитывая (3) и свойство гotgгad ф = 0, получаем:

((^2є + і^а) Е,gгadф) = 0, Vф Є Н0, (10)

((^2є + і^а) Е, gгadф) = J (^2є + і^а) Е ■ gгadф^П,

п

J (^2є + і^а) Е ■ gгadф^П = у div [(^2є + і^а) Е] ф^П. (11)

пп

Из (11) следует, что уравнение (10) является вариационным аналогом закона сохранения (2), т.е. решение вариационной задачи (9) удовлетворяет закону сохранения заряда (2) в слабом смысле.

Из (8) следует, что любая векторная функция Е Є Н0(го!; П) может быть представлена с использованием декомпозиции Гельмгольца Е = и + gгadф, где и Є Н0(го!; П) П

Н0^іу; П), ф Є Н0(П) [8].

3. Дискретные аналоги вариационных задач

Для построения дискретного аналога вариационной задачи элементы пространства Н(го!; П) аппроксимируются элементами дискретного подпространства Н2(го!; П). Введем дискретный аналог вариационной задачи.

Дискретная вариационная задача. Для J0re Є Ь2(П)3 найти Е^є Є Н(го!; П) и Е2т Є Н2(го!; П) такие, что для Є Н^(го!; П) и Є Н^(го!; П):

(і-1У X Е^е, Ух ^)п — (^Е*, ^)п — (иаЕ^, ^)п = 0,

(і 1у х V х ^^П — (w2єEhm, '^П + (^аEhe, ^)п = (WJ0re,

Поскольку для построенных дискретных подпространств имеет место свойство включения

и2 Є Hh(gгаd; П) ^ Уи2 Є Н2(го!; П),

аппроксимация напряженности электрического поля Е2 удовлетворяет закону сохранения заряда в слабой форме:

(—^Е2т — аЕ2е, Уи2) = 0, Є Hh(gгad; П),

(^Е2е — аЕ2т, Уи2) = 0, ^2 Є Н2^; П).

Раскладывая действительную и мнимую компоненты поля Е2 по базису дискретного подпространства Н^(го!; П) и выбирая в качестве функций VI и v2 базисные функции wm, перейдем к эквивалентной системе линейных алгебраических уравнений:

(12)

где Ег и Е — веса в разложении по базису действительной и мнимой компонент поля Е2 соответственно. Элементы матриц В, С, О и вектора ^ определяются следующими соотношениями:

[О]г- =( -^го! w2, го! w2

В + В — с Ег 0

С В + В . Еі . /

.2 2 „т2\

і

[сЬ = (а^2, ^2) п,

[/]і = (^0ге, w2)п• Определяется следующуя матрица предобусловливания:

-1

(В + В )і

Сіі

—Сіі (Я + В )і

(13)

где (Я + В)іі и Сіі — главные диагонали соответствующих матриц.

4. Локальные векторные базисные функции на тетраэдральной сетке

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

Для векторной величины и определим на тетраэдре К интерполяционные функции и2 в виде

и2 (х) = £ «,(иК2(х), (14)

гей1

где — степени свободы, Б — множество индексов степеней свободы, w2 — базисные функции, х = (х1,х2,хз) — координата в Лз.

Пусть рг(г = 1, 2, 3, 4) — координаты вершины произвольного тетраэдра К (рис. 3), Аг(х) — трехмерные барицентрические координаты точки х € К относительно вершин тетраэдра. Тогда локальные векторные базисные функции первого порядка первого типа имеют вид

wf = А1УА2 — А2УАь wK = А^Аз — Аз^А1,

wK = А1УА4 — А4УА1, wK = А2УАз — АзУА2, (15)

wK = А2VА4 — А4^А2, w(f = АзУА4 — А4УАз.

Векторные базисные функции (15) ассоциированы с ребрами тетраэдра. Связь локальной нумерации ребер с нумерацией вершин тетраэдра показана на рис. 3. Пространство базисных функций первого типа обозначим Н2(го!; П; 1).

Р4

Рис. 3. Нумерация вершин и ребер тетраэдрального конечного элемента

Выполнение закона сохранения зарядов в слабой форме (10) связано со свойством включения пространств. Покажем, что оно имеет место для дискретного подпространства Н2(го!; П; —). В качестве скалярных локальных базисных функций дискретного подпространства возьмем барицентрические координаты:

фК = Аь фК = Л2,

фК = ^ фК = А4 •

(16)

Дискретное скалярное подпространство, построенное с использованием локальных базисных функций (16), обозначим Н2^га^ П; 1).

Имеют место следующие соотношения:

УфК = УфК УфК

*К

WK — wK + w

К

УА 1 = — w К — w;K — w3 = УА2 = УАз = УА4

W.K + wK — w

К

WK — wK + w

К

'К, 'К,

К

(17)

Из (17) следует, что градиенты скалярных базисных функций пространства Н2^га^ П; 1) являются линейной комбинацией векторных базисных функций пространства Н2(го!; П; 1). То есть имеет место свойство включения:

и € Н2^га^ П; 1) ^ Уи € Н2(го!; П; 1).

Определим локальные векторные базисные функции первого порядка второго типа, образующие базисные функции пространства Н2(го!; П; 2):

w

w

w

w

w

w

К

1,1

К

2,1

К

3.1 К

4.1 К

5.1 К

6.1

А1УА2,

А1УА3,

А1УА4,

А2УА3,

А4УА2,

А3УА4,

w

w

w

w

w

W

К

1,2

К

2,2

К

3.2 К

4.2 К

5.2 К

6.2

А2 УА1,

А3 VА1,

А4 V А1, А3 УА2, А2 V А4, А4 УА3.

(18)

С каждым г-м ребром тетраэдра ассоциированы две локальные базисные функции

wi,1 и wi,2.

Дискретному подпространству Hh(rot;^;2), натянутому на векторные базисные функции второго типа, соответствует дискретное подпространство Hh(grad; П; 2) со скалярными локальными базисными функциями второго порядка [9]:

ФK = Лі(2Лі —І), ФK — Л2(2Л2 —І),

ФK = Лз(2Лз —І), Ф^Т = Л4(2Л4 —І),

ФK = 4ЛіЛ2, ФK = 4ЛіЛз,

ФK = 4ЛіЛ4, ФГ = 4Л2Лз,

ФK = 4Л2Л4, ФKo = 4Лз Л4.

Построим иерархический базис пространства Hh(rot; П; 2):

w

w

w

w

w

w

K

і,і

K

3.1

K

5.1 K

1.2 K

3.2

K

5.2

Л^Л2 — Л2VЛl, ЛlVЛ4 — Л4VЛl, Л2VЛ4 — Л4VЛ2,

ЛlVЛ2 + Л2VЛl, ЛlVЛ4 + Л4VЛl, Л4VЛ2 + Л4VЛ2,

w

w

w

w

w

w

K

2,1

K

4.1 K

6.1 K 2,2 K

4.2 K

6.2

Лі VЛз — AзVAl, Л2 VЛз — ЛзVЛ2, Лз VЛ4 — Л4VЛз. Лі VAз + ЛзVЛl, Л2 VЛз + ЛзVЛ2, Лз VЛ4 + Л4VЛз.

(19)

Иерархический базис пространства Hh(grad; П; 2) имеет следующий вид:

ФK = Лі, K2 = Л2,

ФK = Лз, Ф4K = Л4,

ФK = ЛіЛ2, ФK = ЛіЛз

ФK = ЛіЛ4, ФГ = Л2Лз

ФK = Л2Л4, Фй = Лз Л4

В дальнейшем в работе рассматриваются иерархические базисы пространств Hh(rot; П; 2) и Hh(grad; П; 2) на тетраэдральной сетке.

5. Двухуровневый итерационный решатель

Рассмотрим систему линейных алгебраических уравнений с невырожденной матрицей А размерности п х п:

Ах = Ъ. (20)

Обозначим через V подпространство пространства И”- с dim(V) = т. Введем матрицу Р размерности п х т, столбцы которой являются базисом подпространства V:

Р IV ^ И”.

Сформулируем двухуровневый итерационный алгоритм решения СЛАУ, использующий подпространство V:

Алгоритм SV(A,b,xov) ro = b - Axo; for i = 1, 2,... g = PT ri-i;

z = (PTAP)-1g or z = S((PTAP),g,Yp);

Xi—1/2 = Xi-i + Pz; ri—1/2 = b — Axi—1/2; y = S (A,ri—1/2,7); xi = xi—1/2 + y; ri = b - Axi;

if ||ri| < v||b||, then return xi; increase i,

где S(A, b, y) — итерационный решатель для СЛАУ Ax = b.

Матрица, полученная после дискретизации рассматриваемой задачи, не обладает симметрией, поэтому в качестве решателя будет применяться предобусловленный стабилизированный метод бисопряженных градиентов BICGSTAB(A, b, y) [10].

Основываясь на работе R. Hiptmair [11], в качестве подпространства V в алгоритме Sv(A, b,x0v) можно выбрать линейное пространство, эквивалентное ядру rot оператора

Nh(rot; П) = {u 6 Hh(rot; П); V х u = 0}.

Столбцы матрицы P — координаты градиентов базисных функций пространства H(grad; П) в базисе пространства H(rot;n). В таком случае система линейных алгебраических уравнений

PT APu = PT r (21)

эквивалентна следующей дискретной вариационной задаче.

Дискретная вариационная задача. Для F0re, Foim 6 Ь2(П)3 найти U^, 6 H(grad; П) и 6 H[j(grad; П) такие, что для VvJ1 6 H[j(grad; П) и Vv^ 6 H(grad; П) выполняются равенства:

(^—1V х VUhe, V х Vvh)n - (c^VUe, Vvh)n - (caVU^, V<)n = (Foim, Vv2h)n,

(^—V х VUim, V х Vvh)n - (c2eVUfm, Vvh)n + (caVU^e, Vv2h)n = (Fore, Vv2h)n.

Учитывая, что V х V = 0, получаем дискретную вариационную задачу.

Для Fore, Foim 6 Ь2(П)3 найти U^ 6 H^(grad; П) и U^ 6 H^(grad; П) такие, что для Vv ^ 6 H^(grad; П) и Vv^ 6 H^(grad; П) выполняются равенства:

-(c2eVUrhe, Vv h)n - (caVUt, Vv *)n = (Foim, Vv2h)n,

-(c2eVUt, Vv2h)n + (caVUh, Vv2h)n = (Fore, Vv2h)n.

Таким образом, матрица системы (21) может быть получена с помощью конечноэлементной технологии, а не путем непосредственного умножения матриц. Решатель SV(A,b,xov), использующий ядро rot, обозначим SN(rot)(A, b, xov).

6. Мультипликативный алгоритм

Пусть существуют два подпространства V1 и V2 пространства Rn такие, что

V1 U V2 = Rn, V1 П V2 = 0.

Введем матрицы P1 и P2, столбцы которых образуют базисы подпространств V1 и V2 соответственно.

Для последовательного уточнения начального приближения элементами подпространств V1 и V2 можно сформулировать итерационный алгоритм решения СЛАУ (20)

[12].

Мультипликативный алгоритм Шварца: ro = b - Axo; для i = 1, 2,... g1 = P1T ri—1;

Z1 = (PTAP1)—1g1 или Z1 = S1((PTAP1),g1 ,Y1); xi—1/2 = xi—1 + P1z1; ri—1/2 = b - Axi—1/2; g2 = P2T ri—1/2;

Z2 = (PT AP2) —1g2 или Z2 = S2 ( (PT AP2) , g2 ,Y2); xi = xi—1/2 + P2Z2; ri = b - Axi;

если xi удовлетворяет необходимой точности, то закончить; увеличить i.

Введем следующие пространства:

Nh(rot; П; 2) = {u 6 Hh(rot; П;2); V х u = 0},

Nh(rot; П; 1) = {u 6 Hh(rot; П; 1); V х u = 0}.

Согласно определению, иерархический базис пространства Hh(rot; П; 2) состоит из базисных функций пространства Hh(rot; П; 1) и градиентов скалярных функций.

Следовательно, Hh(rot; П; 2) можно представить таким образом:

Hh(rot; П; 2) = Hh(rot; П; 1) U Nh(rot; П; 2),

Hh(rot; П; 1) П Nh(rot; П; 2) = Nh(rot; П; 1).

Пусть V1 = H^rot;^!), V2 = Nh(rot;П;2), тогда мультипликативный алгоритм Шварца можно использовать для решения СЛАУ (12) [13].

Из определения подпространств следует, что матрицы (P^AP1) и (P2TAP2) будут соответствовать дискретным вариационным задачам, сформулированным относительно векторного уравнения Гельмгольца в пространствах Hh(rot; П; 1) и Nh(rot; П; 2) соответственно. Таким образом, для получения данных матриц будем использовать конечноэлементную технологию.

В качестве решателя S1((P1TAP1),g1,Y1) используем стабилизированный метод би-сопряженных градиентов. Согласно определению иерархического базиса пространства Hh(grad; П; 2), имеет место следующее свойство:

Nh(rot; П; 1) С Nh(rot; П; 2).

Таким образом, решатель S2((P2tAP2), g2, y2) строится на основе двухуровневого алгоритма SV(A,b, 0,y1 ), где V = Nh(rot; П; 1).

7. Результаты численного тестирования алгоритма на модельной задаче

Численное исследование предложенного алгоритма выполнялось на модельной задаче

(1)-(6) моделирования гармонического электрического поля в области [—1,1]3, состоящей из двух подобластей с различными значениями удельной проводимости и а2 (рис. 4). В подобласти 1 расположена генераторная катушка, в которой задан ток силой 1 А с частотой 14 МГц. Границей между подобластями является плоскость х = 0.

Таблица 1. Длина наибольшего ребра сетки Ътах и размерности дискретных подпространств, построенных на последовательности неструктурированных тетраэдральных сеток

Сетки ъ ^тах И(го^ 0; 1) Х(го^ 1) И(гс^; 0; 2) Х(го^ 2)

Т1 2 5 1 О 1 8776 1474 17552 10250

Т2 8.84 ■ 10-2 63756 9632 127512 73388

Тз 6.25 ■ 10-2 236364 34204 472728 270568

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

Сетки (^1> 02)

(1,1) (0,0) (1,10) (1,0.1) (1,0) (10,1) (0.1,1) (0,1)

Т1 4 3 5 5 5 4 5 6

Т2 64 53 92 93 118 51 93 119

Т3 767 544 1101 1457 1879 590 1545 1702

Таблица 3. Время, с, нахождения приближенного решения алгоритмом ^(гсЛ)^, Ь, Жо^) при различных соотношениях проводимости на различных тетраэдральных сетках

Сетки (01> 02)

(1,1) (0,0) (1,10) (1,0.1) (1,0) (10,1) (0.1,1) (0,1)

Т1 6 3 3 8 8 5 8 10

Т2 122 56 118 196 215 71 205 269

Т3 2651 1005 1674 6957 7224 1364 5774 7022

Тестирование проводилось на компьютере с процессором Athlon-XP +1800 и 512 Мб оперативной памяти.

В табл. 1 указаны размерности дискретных подпространств, построенных на последовательности неструктурированных тетраэдральных сеток Т, а также длина наибольшего ребра сетки Нтах.

Время (в секундах), необходимое мультипликативному алгоритму и алгоритму ^(го-ь) (А, Ь, х0V) при уменьшения относительной нормы невязки в 107 раз, показано для различных соотношенияй между значениями удельной проводимости в первой и второй подобластях на различных тетраэдральных сетках в табл. 2 и 3 соответственно.

Серия векторных суперкомпьютеров NEC SX началась в 1983 году с модели SX-1/2, первого суперкомпьютера, предоставлявшего пиковую производительность, превосходящую 1 ГФлопс. Модель NEC SX-8 стала доступной в январе 2005 года, а система, установленная в HLRS, начала работать в апреле 2005 года Подробное описание системы находится в [14]. Код для математического моделирования прибора ВИКИЗ (далее просто код ВИКИЗ) был изначально написан для работы на РС. Поэтому требовалось портирование кода на NEC SX-8 (рис. 5).

Сетка строилась с помощью автоматического З^генератора сеток NETGEN [15]. К сожалению, из-за проблем при компиляции библиотек Тс1, Тк и Тіх нам не удалось скомпилировать NETGEN для запуска на SX-8, поэтому NETGEN был скомпилирован для запуска на фронтальном узле. На данном этапе это не создало проблем, поскольку сетка строилась до начала работы решателя и не изменялась в процессе его работы. Тем не менее, распараллеливание процесса построения сетки остается задачей на будущее, так как время работы последовательной версии составляет около 5 ч.

8. Расчет ВИКИЗ на N£0 ЯХ-8

г

Рис. 5. Узел №0 SX-8

Рис. 6. Зависимость разности фаз от глубины зондирования

При портировании решателя никаких проблем не возникло.

Первые тесты последовательной версии кода ВИКИЗ были проведены для проверки корректности портирования кода на SX-8 с положительным результатом.

На втором шаге восемь копий последовательной программы одновременно запускались на одном узле SX-8. После отработки всех восьми копий запускались следующие восемь копий и т. д. Такой подход вполне рационален, поскольку входные данные состоят из независимых файлов.

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

На рис. 6 показана численно полученная зависимость разности фаз от глубины зондирования: z — глубина зондирования; Аф — разность фаз, которая позволяет геофизикам определить свойства породы. Эта величина результат реальных измерений зондом ВИКИЗ.

9. Перспективы

Следующие вопросы будут рассмотрены в процессе дальнейшей работы над численным моделированием ВИКИЗ:

— переход от NETGEN к другому генератору сеток для параллельного построения сетки;

— переход от внешнего параллелизма к внутреннему, поскольку алгоритм позволяет эффективное распараллеливание;

— переход от гармонических полей к полям, меняющимся со временем;

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

Список литературы

[1] Технология исследования нефтегазовых скважин на основе ВИКИЗ: Методическое руководство / Ред. М.И. Эпов, Ю.Н. Антонов. Новосибирск: НИЦ ОИГГМ СО РАН, Изд-во СО РАН, 2000.

[2] Еремин В.Н. Прибор высокочастотного электромагнитного каротажа в процессе бурения ВИКИПБ-7 // Геофизический вестник. 2005. № 1. C. 15-19.

[3] HlPTMAlR R. Finite elements in computational electromagnetism // Acta Numerica. 2002. Vol. 11. P. 237-339.

[4] Nedelec J.-C. Mixed finite elements in R3 // Numerische Mathematik. 1980. Vol. 35. P. 315-341.

[5] Nedelec J.-C. A new family of mixed finite elements in R3 // Numerische Mathematik. 1986. Vol. 50. P. 57-81.

[6] Monk P. Finite element methods for Maxwells equations. Oxford: Clarendon Press, 2003.

[7] Arnold D.N., Falk R.S., Winther R. Finite element exterior calculus, homological techniques an applications // Acta Numerica. 2006. P. 1-155.

[8] Caorsi S., Fernandes P., Raffetto M. On the converges of Galerkin finite element approximation of electromagnetic eigenproblem // SIAM J. Numer. Anal. 2000. Vol. 38. P. 580-607.

[9] Webb J.P. Hierarchal vector basis functions of arbitrary order for triangular and tetrahedral finite elements // IEEE Trans. Antennas Propag. 1999. Vol. 47, N 8. P. 1244-1253.

[10] Van der Vorst H.A. BiCGSTAB a fast and smoothly converging variant of BiCG for the solution of nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1992. Vol. 13, N 2. P. 631-644.

[11] Hiptmair R. Multigrid methods for Maxwell’s equations // SIAM J. Nymer. Anal. 1998. Vol. 36, N 1. P. 204-225.

[12] Saad Y. Iterative methods for sparse linear systems. PWS Publishing Company, 1996.

[13] Nechaev O.V., Shurina E.P., Botchev M.A. Multilevel iterative solvers for the edge finite element solution of the 3D Maxwell equation // Department of Applied Mathematics Internal Report 1806, University of Twente, 2006.

[14] Tagaya S., Nishida M., Hagiwara T. et al. The NEC SX-8 Vector Supercomputer System / Ed. M. Resch, T. Bonisch, K. Benkert, T. Furui, Y. Seo, W. Bez. High Performance Computing on Vector Systems 2005. Proceedings of the High Performance Computing Center Stuttgart, March 2005. Springer, 2006. P. 3-24.

[15] http://www.hpfem.jku.at/netgen/

Поступила в редакцию 19 ноября 2007 г.