ФИЗИКА ЗЕМЛИ, 2014, № 4, с. 50-55
УДК 550.831:517.968:519.612:004.021
ИЗУЧЕНИЕ СТРУКТУРНЫХ ОСОБЕННОСТЕЙ ГРАВИТАЦИОННОГО И МАГНИТНОГО ПОЛЕЙ ЛИТОСФЕРЫ С ИСПОЛЬЗОВАНИЕМ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ
© 2014 г. П. С. Мартышко1,3, Н. В. Федорова1, Е. Н. Акимова2,3, Д. В. Гемайдинов1,3
Институт геофизики им Ю.П. Булашевича УрО РАН, г. Екатеринбург E-mail: mpdir@mail.ru 2Институт математики и механики УрО РАН, г. Екатеринбург 3Уральский федеральный университет им. Б.Н. Ельцина, г. Екатеринбург E-mail:pmart3@mail.ru, nataliavf50@mail.ru, aen15@yandex.ru, gemaidinov@gmail.com
Поступила в редакцию 26.03.2013 г.
В статье описываются параллельные алгоритмы для изучения структурных особенностей аномалий гравитационного и магнитного полей литосферы, основанные на повысотных трансформациях данных. Алгоритмы численно реализованы на суперкомпьютере "Уран". С помощью разработанной компьютерной технологии построены карты региональных и локальных аномалий гравитационного и магнитного полей для северо-восточного сектора Европы в пределах трапеции 48°—62° в.д. и 60°—68° с.ш.
DOI: 10.7868/S0002333714040097
ВВЕДЕНИЕ
В геофизической практике для разделения потенциальных полей на составляющие широко используются методы пересчета в верхнее и нижнее полупространство. Приведение полученных результатов на начальный уровень дает возможность разделять длинноволновые и коротковолновые составляющие амплитудного спектра аномалий, созданных литосферой.
В работах [Мартышко, Пруткин, 2003; Мартышко и др., 2010] для разделения источников гравитационного поля по глубине была предложена следующая методика трансформации полей. На начальном этапе решается задача о выделении эффекта источников в слое от земной поверхности до некоторой глубины Н. Поле продолжается вверх на уровень Н, при этом влияние локальных приповерхностных источников (до глубины Н) если и не устраняется совсем, то значительно ослабевает. Для того чтобы окончательно избавиться от влияния локальных источников, находящихся в горизонтальном слое от дневной поверхности до глубины Н, пересчитанное вверх поле затем продолжается вниз на глубину Н. При этом, поскольку задача относится к классу некорректно поставленных задач необходимо использовать методы с применением регуляризации. На последнем шаге поле пересчиты-вается вновь вверх на уровень дневной поверхности. Полученное поле можно рассматривать как поле от источников, расположенных ниже границы Н. Далее вычитаем это поле из наблюденного
и получаем поле от слоя. Повторяя эту процедуру для различных значений высот и глубин, выделяем поля слоев с соответствующими границами.
Такая методика разделения гравитационного поля была использована при построении плот-ностной модели земной коры Уральского региона [Мартышко и др., 2010; 2011]. При изучении крупных территорий приходится задавать большие массивы данных, что приводит к значительным затратам времени при вычислениях на однопроцессорных компьютерах. Использование параллельных алгоритмов для многопроцессорных вычислительных систем значительно сокращает время расчетов.
В данной работе приведено описание математического аппарата и алгоритмов параллельных вычислений, использованных при создании компьютерной технологии и результатов применения этой технологии для изучения структурных особенностей гравитационного и магнитного полей литосферы северо-восточного сектора Европы.
МАТЕМАТИЧЕСКИЙ АППАРАТ, ИСПОЛЬЗУЕМЫЙ ДЛЯ ТРАНСФОРМАЦИЙ
Введем декартову прямоугольную систему координат с осью Z, направленной вниз, плоскость ХОУ совпадает с дневной поверхностью. Пересчет измеренного на площади Б =
= {(х, у) е Я2:а < х < Ь, с < у < й} земной поверхности поля и (х, у, г) г=о вверх на уровень г = — Н реализуется по формуле Пуассона
Ux у, -H) = 2П ii
H
2п-1-1 {(х - X')2 + (у - у) + и2]3/2 (1) х и(X, y',0)dx,dy,. Для пересчета поля вниз на глубину г = Н и нахождения значений и(х, у, -И) решается уравнение Фредгольма первого рода
Ku =
i ir
H
2п -1-1 [(x - x')2 + (y - у')2 + H2 ]
х u(x', y',0)dx'dy' = U(x, y, H)
3/2
(2)
k+1 k U = u -
(K (Kuk - U), Kuk - U)
IK(Kuk - U)
3) метод минимальной ошибки
(Kuk - U); (5)
k+1 k u = u --
Kuk - U
\KT (uk - U) 4) метод наискорейшего спуска
KT (uk - U); (6)
k+1 k u = u -■
KT (Ku - U)
KKT (uk - U)
KT (uk - U). (7)
Заметим, что для методов (5)—(7) в регуляри-зованном варианте матрица К заменяется на К + а I.
Условием обрыва итерационных процессов (4)— (7) является выполнение условия
\\Ku - U
< 6
Уравнение 1-го рода (2) представляет собой некорректную задачу, поэтому необходимо при его решении использовать методы регуляризации. Оператор уравнения (2) является положительно определенным и самосопряженным, поэтому можно применять схему М.М. Лаврентьева [1962].
После дискретизации на сетке уравнения (2) и аппроксимации интегрального оператора по квадратурным формулам задача сводится к решению системы линейных алгебраических уравнений (СЛАУ) с симметричной матрицей К. В регу-ляризованном виде система имеет вид
(К + а 1)и = и, (3)
где I — единичный оператор, а — параметр регуляризации.
Для решения СЛАУ могут быть использованы итерационные методы градиентного типа [Васин, Еремин, 2005]:
1) итеративно регуляризованный метод простой итерации
ик+1 = ик--1— [(К + а1) ик - и], (4)
^ тах
где ^тах — максимальное собственное значение матрицы К + аI (симметричный случай);
2) метод минимальных невязок
при достаточно малом s.
Эти методы успешно применяются на многопроцессорных вычислительных системах различного типа при решении обратных задач гравиметрии и магнитометрии для восстановлении поверхности раздела между слоями и плотности в слое [Акимова, Гемайдинов, 2007; Акимова, 2009; Мартышко и др., 2010].
ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ И ОЦЕНКА МАСШТАБИРУЕМОСТИ
Параллельные алгоритмы повысотных трансформаций на основе итерационных методов градиентного типа разработаны и численно реализованы на суперкомпьютере "Уран", установленном в Институте математике и механике УрО РАН.
Суперкомпьютер "Уран" — многопроцессорный вычислитель кластерного типа, каждый вычислительный узел которого состоит из двух 4-ядер-ных процессоров Intel Quad-Core Xeon, работающих на частоте 3.00 ГГц, и содержащих 16/32 Гбайт оперативной памяти. В общей сложности пользователю доступны 1664 вычислительных ядра и 3584 Гбайт оперативной памяти.
Распараллеливание итерационных методов градиентного типа (4)—(7) (этап пересчета поля вниз) основано на разбиении матрицы K горизонтальными или вертикальными полосами на m блоков, а вектора решения u и вектора правой части U СЛАУ на m частей так, что n = m х L, где n — размерность системы уравнений, m — число процессоров, L — число строк или столбцов матрицы в блоке.
На текущей итерации каждый из m процессоров вычисляет свою часть вектора решения. В случае умножения матрицы K на вектор u каждый из m процессоров умножает свою часть строк матрицы K на вектор u. Host-процессор (ведущий) отвечает за пересылки данных и также вычисляет свою часть вектора решения.
На этапе пересчета поля вверх при вычислении интегрального оператора (1) используется параллельный алгоритм умножения матрицы на вектор.
Программы, реализующие параллельные алгоритмы повысотных трансформаций на суперкомпьютере "Уран", написаны с помощью библиотеки MVAPICH2 (MPI) на языке программирования Фортран.
х
2
52
МАРТЫШКО и др.
Таблица
m (число процессоров) T m (время, мин) Mm (масштабируемость)
50 3.73 -
100 1.88 1.98
150 1.27 2.95
300 0.67 5.62
Заметим, что для решения задачи с двойной точностью итерационными методами на суперкомпьютере "Уран" для хранения матрицы K требуется достаточно много памяти. Для сетки 300 х 300 количество элементов матрицы K занимает примерно 60.3 Гб, а каждому отдельному процессу MPI — программы системой запуска задач на кластере "Уран" выделяется стандартно не более 2 Гб. Поэтому предварительно матрица K разбивается вертикальными полосами на блоки и хранится в памяти каждого процессора по частям. Таким образом, с учетом требуемой памяти для решения нашей задачи необходимо не менее 50 процессоров "Уран".
Для сравнения времени решения задачи введем коэффициенты ускорения и эффективности параллельных алгоритмов
Sm = TjTm, Em = sjiï, S = TJT2, (8)
где Tm — время выполнения параллельного алгоритма на суперкомпьютере "Уран" с числом процессоров m (m > 1), T1 — время выполнения последовательного алгоритма на одном процессоре. Tm представляет собой совокупность чистого времени счета и накладных расходов на межпроцессорные обмены Tm = Tc + To.
В связи с ограничением по памяти эффективное решение задачи на одном процессоре "Уран" невозможно. При использовании технологии записи данных на жесткий диск и скачивания с жесткого диска результирующего вектора решения uk на каждой итерации МПИ время счета увеличивается более чем в 1000 раз и составляет для сетки 300 х 300 более 62 ч.
Поскольку в нашем случае оценить эффективность решения задачи по формуле (8) нельзя, то в качестве показателя эффективности предлагается использовать понятие масштабируемости параллельного алгоритма.
Под коэффициентами масштабируемости параллельного алгоритма мы будем понимать возможность ускорения вычислений в M раз пропорционально числу используемых процессоров m с учетом того, что объем данных не превосходит ограничений по памяти для каждого процессора
Мт Тшт! Тт, (9)
где Ттп — время выполнения параллельного алгоритма на минимально возможном числе процессоров с учетом объема выделяемой памяти.
В таблице приводится время выполнения одной повысотной трансформации поля (решение последовательно задач (1), (2) и снова (1) в соответствии с методикой) на суперкомпьютере "Уран" с помощью параллельных алгоритмов ПАМВ (параллельный алгоритм умножения матрицы на вектор) и ПМПИ (параллельный метод простой итерации) и значение коэффициентов масштабируемости в зависимости от количества использованных процессоров т.
Вычислительные эксперименты показали, что параллельные алгоритмы имеют высокую степень масштабируемости, поскольку с увеличением числа процессоров в М раз время счета уменьш
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.