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

Аннотация
1. Линейная обратная задача гравиметрии о нахождении плотности в слое.
Для решения линейной обратной задачи гравиметрии предложены и численно
реализованы на суперкомпьютере МВС─1000 и видеоускорителях NVIDIA GeForce
Параллельные итерационные алгоритмы (МПИ и ММН).
Решена обратная задача гравиметрии с реальными данными.
2. Задачи электроразведки (СЛАУ с блочно-трехдиагональными матрицами).
Для решения СЛАУ предложены и численно реализованы на видеоускорителе NVIDIA
GeForce, 4-х ядерном процессоре Intel Core I5-750 и суперкомпьютере МВС-1000
параллельный алгоритм матричной прогонки и параллельный метод сопряженных
градиентов с предобуславливателем. Решена модельная задача о нахождении
распределения потенциала на поверхности земли в проводящей среде.
Проведено сравнение времени счета параллельных алгоритмов на вычислителях
различного типа с анализом эффективности и ускорения.
Работа выполнена в рамках Междисциплинарного проекта УрО РАН и Программы
фундаментальных исследований Президиума РАН № 14 при поддержке УрО РАН.

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

Нахождение плотности в слое
Задача сводится к решению линейного двумерного интегрального уравнения Фредгольма первого рода для нахождения искомой плотности [1]
(1)
где гравитационная постоянная, гравитационный эффект,
порождаемый источниками в горизонтальном или криволинейном слое.
Уравнение гравиметрии (1) относится к классу некорректно поставленных задач, решение которой обладает сильной чувствительностью к погрешности правой части, полученной в результате измерений и предварительной обработки геофизических данных.
____________________________________________________________________
[1]. Martyshko P.S., Koksharov D.E. On the construction of the density sections using gravity data
// Extended Abstracts of 66th EAGE Conference and Exhibition. Paris, 7-12 June 2004. P. 143.

Методы решения
После дискретизации уравнения (1) на сетке где задана правая
часть и аппроксимации интегрального оператора по квадратурным
формулам задача (1) сводится к решению СЛАУ с плохо обусловленной
либо симметричной матрицей (горизонтальный слой), либо несимметричной
матрицей (криволинейный слой)
В случае криволинейного слоя СЛАУ предварительно преобразуется к виду
где параметры регуляризации.
Для решения СЛАУ (2), (3) используются регулярные итерационные методы градиентного типа [2].
______________________________________________________________________
[2]. Васин В.В., Еремин И.И. Операторы и итерационные процессы Фейеровского типа.
Теория и приложения. – Екатеринбург, 2005.

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

Параллельная реализация на МВС-1000
Ранее для решения линейной обратной задачи гравиметрии о восстановлении
переменной плотности в слое параллельные методы градиентного типа были численно
реализованы на МВС-1000 с высокой эффективностью распараллеливания [3].
Многопроцессорный вычислительный комплекс МВС─1000 ─ российский массивно-
параллельный суперкомпьютер кластерного типа с распределенной памятью,
установленный в ИММ УрО РАН.
МВС1000/64 (UM64): 14 2-х процессорных 2-х ядерных модулей AMD
Opteron 64 bit (2.6 ГГц); интерфейс GbitEthernet; 112 Гб оперативной
памяти.
______________________________________________________________________________
[3]. Акимова Е.Н., Гемайдинов Д.В. Параллельные алгоритмы решения задачи гравиметрии
о восстановлении плотности в слое // Труды института математики и механики УрО РАН. –
Екатеринбург: УрО РАН, 2007. Т. 13. № 3. С. 3-21.

Рис. 2. Схема распределения данных по процессорам
Комплекс параллельных численных алгоритмов реализован с помощью библиотеки MPI на языке Фортран [4]. Распараллеливание итерационных методов основано на разбиении матриц и горизонтальными полосами на блоков, а вектора решения и вектора правой части СЛАУ на частей так, что где размерность системы уравнений,
число процессоров, число строк матрицы в блоке.
Параллельная реализация на МВС-1000
______________________________________________________________________
[4]. Баранов А.В., Лацис А.О., Сажин C.В., Храмцов М.Ю. Руководство пользователя системы МВС-1000.
URL: http://parallel.ru/mvs/user.html.

На МВС-1000/64 решена задача с реальными данными о восстановлении плотности
в горизонтальном слое между глубинами для области
Шаги сетки: Гравитационная постоянная
Задача сводится к СЛАУ с плохо обусловленной симм. заполненной матрицей
Для решения задачи использовался параллельный итеративно регуляризованный МПИ
с параметром регуляризации
Рис. 3. Исходное гравитационное поле для области
Решение задачи гравиметрии с реальными данными

Рис. 4. Линии уровня и распределение восстановленной
плотности в слое 10 - 20 км для области
Более темные участки соответствуют зонам разуплотнения,
представляющим интерес для геофизической интерпретации.

В табл. 1 приведены времена счета на МВС-1000/64 и коэффициенты ускорения
и эффективности решения задачи гравиметрии для области с использованием параллельного МПИ для 200  200 точек сетки (430 итер).
Матрица СЛАУ формируется и хранится в памяти каждого процессора по частям.
Таблица 1. Решение задачи гравиметрии о восстановлении плотности в слое
Ускорение и эфективность
где время выполнения последовательного алгоритма на одном процессоре,
время выполнения параллельного алгоритма на МВС с числом процессоров
совокупность чистого времени счета и накладных расходов.

Распараллеливание на видеоускорителях с помощью технологии CUDA http://www.nvidia.ru
Для организации параллельных вычислений актуальным в настоящее время
является использование видеоускорителей (GPU) компании NVIDIA (рис.4).
Базовый блок - мультипроцессор, содержащий процессорные ядра.
Работа нескольких ядер мультипроцессора основана на архитектуре типа SIMD.
Для поддержки параллельных вычислений компания NVIDIA разработала технологию CUDA [5] - среду разработки программ на языке Cи, позволяющую создавать программное обеспечение для решения сложных вычислительных задач.
Линейная задача гравиметрии решена на видеоускорителях GeForce GTX 285 (GPU-1)
и GeForce GTX 260 (GPU-2) с помощью технологии CUDA .
Рис. 5. Видеоускоритель GeForсe GTX 285
____________________________________________________________
[5]. Берилло А. NVIDIA CUDA – неграфические вычисления на графических процессорах.
URL: http://www.ixbt.com/video3/cuda-1.shtml

Табл. 2. Техн. характеристики вычислительной платформы

В табл. 3 приводятся времена решения задачи гравиметрии на Intel Core I5-750 без использования и с использованием видеоускорителей GeForce на разных сетках: 110 x 110 (матр. СЛАУ 12100 x 12100) и 200 x 200 (матр. СЛАУ 40000 x 40000).
Результаты численных экспериментов
Сравнение времени счета на GPU-1,2 и МВС-1000/64

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

1. Задача ВЭЗ о нахождении потенциала на поверхности земли сводится
к решению двумерного уравнения Лапласа в цилиндрической системе координат [7]
(6)
с граничными условиями непрерывности потенциала и непрерывности нормальной
составляющей к границам плотности тока
условия на горизонтальных прямых (7)
слева,
сверху,
справа и снизу.
После использования конечно-разностной аппроксимации краевая задача (6)–(7) сводится к решению СЛАУ
с блочно-трехдиагональной матрицей.
2. Задача бокового каротажного зондирования (БКЗ) – при измерении потенциала электрического поля использует однотипные зонды разной длины. В результате интерпретации данных каротажа получают значение удельного электрического сопротивления пласта, близкое к истинному, по величинам которого выявляют полезные ископаемые. После конечно-разностной аппроксимации задача БКЗ сводится к решению СЛАУ
с блочно-трехдиагональной матрицей [8]
матрица размерности с блоками Рис.7.
_____________________________________________________________________
[7]. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1966.
[8]. Дашевский Ю.А., Суродина И.В., Эпов М.И. Квазитрехмерное мат. моделирование диаграмм неосесимм.
зондов постоянного тока в анизотропных разрезах // Cиб. ЖИМ. 2002. Т. 5. №3 (11). С. 76-91.
коэф. электропров.

где и искомые и заданные векторы размерности n, квадратные матрицы порядка n.
Будем предполагать, что исходная область P – прямоугольник. Разобьем P на L подобластей вертикальными линиями так, что В качестве параметров выберем векторы связывающие неизвестные на сетке по вертикали. В подобластях, определяемых интервалами (K,K+M), рассмотрим задачи
(10)
(11)
Рис. 8.
(12)
_____________________________________________________________________
[9]. Акимова Е.Н. Распараллеливание алгоритма матричной прогонки //
Математическое моделирование. М.: Наука, 1994. Т. 6. № 9. C. 61- 67.
Методы решения СЛАУ с блочно-трехдиагональными матрицами
1. Параллельный алгоритм матричной прогонки
(9)

где и квадратные матрицы порядка n. Схема параллельного алгоритма матричной прогонки: (10)-(11)-(12) → (14) → (13).
Задача (14) решается классическим алгоритмом матричной прогонки.
Задачи (10)-(11)-(12) и (13) решаются независимо на L процессорах.
Утверждение 2. Если для исходной системы (9) выполняются достаточные условия устойчивости
метода матричной прогонки по А.А. Самарскому [10]
причем хотя бы одно из неравенств – строгое, то эти же условия достаточны и для устойчивости метода матричной прогонки при решении системы уравнений (14) относительно параметров [9].
____________________________________________________________________________________________
[9]. Акимова Е.Н. Распараллеливание алгоритма матричной прогонки // Математическое моделирование.
М.: Наука, 1994. Т. 6. № 9. C. 61- 67.
[10]. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.
Утверждение 1. Если решения задач (10), решения задач (11),
решения задачи (12), решения исходной задачи (9) на тогда
После подстановки (13) в (9) получим систему относительно параметров

Для СЛАУ (15) МСГ с предобуславливателем имеет вид [12]
2. Параллельный метод сопряженных градиентов с предобуславливателем
МСГ – быстрый итерационный метод решения СЛАУ с симметричной положительно-определенной матрицей [11].
Введение предобуславливания применяется с целью сущест. ускорения сходимости итерационного процесса.
Исходная СЛАУ заменяется на
Условием выбора предобуславливателя является следующее:
__________________________________________________________________________________________________________
[11]. Фаддеев В.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.: Гос. издат. физ.-мат. литературы, 1963.
[12]. Амосов А.А. Циркулянтно предобусловленный метод сопряженных градиентов и его применение для численного
решения интегрального уравнения переноса излучения // Труды XI Всерос. школы-семинара «Современные
проблемы математического моделирования». Ростов-на-Дону: Издат. Ростов. госун-та, 2005. Выпуск 4. С. 49-65.
условие останова.
Распаралеливание МСГ с предобуславливателем показано на рис. 1 (слайд 8).

Параллельный алгоритм матричной прогонки (ПАМП) и параллельный метод сопряженных градиентов с предобуславливателем (ПМСГ) реализованы на следующих многопроцессорных системах: МВС─1000/64 с помощью технологии MPI, 4-х ядерном процессоре Intel Core I5-750 (CPU)
на языке С++ в среде разработки «Visual Studiо» и видеоускорителе NVIDIA GeForce GTX 285 (GPU)
с помощью технологии CUDA.
Тех. характеристики вычислительных систем приводятся на слайдах 7 и 13.
С помощью методов ПАМП и ПМСГ решена модельная задача о нахождении распределения потенциала в проводящей среде с известным решением (рис. 9).
Исходные данные и модельное решение задачи предоставлены лабораторией скважинной геофизики ИНГГ СО РАН (г. Новосибирск).
Результаты численных экспериментов решения
задачи о нахождении распределения потенциала
Рис. 9.

После дискретизации задача сводится к СЛАУ с плохо обусловленной симметричной положительно-определенной блочно-трехдиагональной матрицей размерности c квадратными блоками порядка 248 (рис. 7, слайд 16).
Приближенное решение задачи сравнивалось с модельным решением с помощью вычисления относительной погрешности
Предварительно находилось число обусловленности матрицы
В случае решения задачи ПМСГ находилось число обусловленности матрицы

Решение задачи методом ПМСГ t=10 час.
Решение задачи методом ПАМП

Основные результаты работы
1. Для решения линейной обратной задачи гравиметрии о восстановлении плотности
в слое предложены и численно реализованы на МВС-1000 и видеоускорителях
NVIDIA GeForce параллельные итерационные алгоритмы. Решена обратная задача
гравиметрии с реальными данными. Проведено сравнение времени
счета параллельных алгоритмов с анализом эффективности и ускорения.
2. Для решения СЛАУ с блочно-трехдиагональными матрицами предложены
и численно реализованы на видеоускорителе NVIDIA GeForce, 4-х ядерном процессоре
Intel Core I5-750 и МВС-1000 параллельный алгоритм матричной прогонки и параллельный метод сопряженных градиентов с предобуславливателем.
Проведено сравнение времени счета параллельных алгоритмов при решении задачи
о нахождении распределения потенциала на поверхности земли в проводящей среде.