CC BY
45
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2013, № 2 (1), с. 160-167
УДК 519.6; 517.958:5
К ВОПРОСУ О МАСШТАБИРУЕМОСТИ ПРОГРАММЫ ДЛЯ ОБРАТНОЙ ЗАДАЧИ ВОЛНОВОЙ ТОМОГРАФИИ
© 2013 г. С.Ю. Романов
Московский госуниверситет им. М.В. Ломоносова
romanov60@gmail .сот
Поступило в редокцию 19.12.2012
Разрабатываются алгоритмы и масштабируемое программное обеспечение на петафлопсных суперкомпьютерах для решения нелинейных коэффициентных обратных задач волновой томографии со сверхбольшим объемом данных и сверхвысоким разрешением. Предложена структура программы, позволяющая распараллеливать вычисления с высокой степенью масштабируемости. Рассмотрены вопросы вычислительной эффективности. Приведены результаты расчетов на 20480 вычислительных ядрах Суперкомпьютерного вычислительного комплекса МГУ «Ломоносов».
Ключевые слово: коэффициентные обратные задачи, волновое уравнение, ультразвуковая томография, параллельные вычисления, масштабируемость, суперкомпьютер.
Введение
В современной медицине чрезвычайно популярны неинвазивные методы диагностики заболеваний, в современной промышленности -неразрушающие методы контроля изделий. Для этих целей интенсивно разрабатываются томографические методы исследования. Высокое разрешение X-ray и MRT томографов обусловлено простыми, адекватными реальности математическими моделями, высокой точностью входных данных. Трехмерные задачи интерпретации данных сводятся в этом случае к решению набора независимых двумерных линейных задач, которые с успехом решаются на персональной ЭВМ. Значительно более скромные результаты в волновой томографии, использующей ультразвуковые, сейсмические, электромагнитные источники зондирования, связаны с тем, что даже в простейших моделях эти задачи являются сложными нелинейными трехмерными задачами.
Известно большое число работ в нашей стране и за рубежом, которые связаны с рассматриваемыми обратными задачами диагностики с использованием волновых источников излучения. Большинство публикаций посвящено решению задач диагностики в линейном или линеаризованном приближениях. Разработано большое количество пакетов программ для акустического и сейсмического зондирования: Schlumberger; Sierra Geophysics (США), пакет программ Quiklog; Merlin Profilers (Research) Ltd (Англия), пакет программ Seismic Kernel System; Professional Geophysics, Inc. (США), пакет программ Refstat. В ультразвуковой
медицине - разработки: TechniScan Medical Systems (США), Karmanos Cancer Institute и Delphinus Medical Technologies, Inc., (США), Institute for Data Processing and Electronics (Г ермания).
В этих разработках обратные задачи, как правило, решаются в лучевом или линеаризованном волновом приближениях [1-3]. К сожалению, линейные модели не всегда хорошо согласуются с нелинейными физическими процессами.
Задачи волновой томографии в рамках нелинейных моделей получили развитие лишь в последние несколько лет. Это связано с высокой сложностью рассматриваемых задач как с математической, так и с вычислительной точки зрения [4, 5]. Существует интегральный подход на основе метода функции Грина для решения обратных задач волновой томографии. В этом подходе для некоторых постановок установлена единственность решения [6], анализируется возможность его численного нахождения с использованием итерационных процедур [7-9].
В предлагаемой работе обратная задача, в отличие от интегрального подхода, рассмотрена как коэффициентная обратная задача для дифференциального уравнения гиперболического типа. Похожие постановки рассмотрены в [10,11].
В настоящей статье рассмотрены некоторые вопросы, относящиеся к разработке программного обеспечения для решения обратных задач волновой томографии. Одна из проблем рассматриваемой задачи состоит в необходимости выполнения огромного объема вычислений. Разработанная программа позволяет решать
обратные задачи волновой томографии в нелинейной постановке на сетках 1000x1000 элементов в томографическом слое. Решение столь масштабных задач приводит к необходимости проведения расчетов на супер-ЭВМ петафлопс-ного уровня производительности. Одна из серьезных проблем в использования супер-ЭВМ, которая будет рассмотрена в работе, состоит в выборе оптимальной схемы распараллеливания задачи для обеспечения сверхмасштабируемости программного обеспечения, повышения эффективности. Полученные результаты являются новыми и находятся на мировом уровне [12]. В работе [13] рассмотрены вопросы оптимизации кода при параллельных расчетах.
Расчеты проводились на кластерной вычислительной системе МГУ «Ломоносов».
Постановка и методы вычислений обратной задачи
Рассмотрим задачу Коши для гиперболического уравнения, которое в скалярном приближении описывает акустическое поле u(r, t) в
Rn x (0,да) (n=2,3) с точечным источником,
располагающимся в точке r0:
c(rК (r,t) - Au(r, t) = 5(r - r0) f (t) , (1)
u(r, t = 0) = ut (r, t = 0) = 0.
Здесь c 05 (r) является скоростью волны в
среде, r £ Rn (n=2,3) - положение точки в пространстве, А - оператор Лапласа по переменной r. Генерируемый источником импульс описывается функцией f (t). Будем предполагать, что неоднородность среды, вызванная только изменениями скорости, располагается в области QcRn (n=2,3), ограниченной поверхностью S, а вне области неоднородности c(r) = c0 = const, где c0 - известна.
Обратная задача состоит в нахождении функции c(r) , r £ Q, описывающей неоднородность, по экспериментальным данным измерения волны U(r, t) на границе S области за
время (0,T) при различных положениях r0 источника для достаточно большого Т.
Обратная задача ставится как задача минимизации квадратичного функционала
и м2
Ф(м(с),с) = ||м \ST -U|| . (2)
Подробно математические аспекты проблемы рассмотрены в работах [11, 14].
В соответствии с классическими традициями томографических исследований обратная задача диагностики 3D объекта рассматривается как
набор двумерных задач. Иная постановка задачи в двумерном случае рассмотрена в работе [15]. В работах [16, 17] решается обратная коэффициентная задача для более сложной модели плоского тела с неизвестными функциями модуля упругости и плотности.
Расчетная модель, на основе которой проведены вычисления, получена из дифференциального уравнения (1). Для решения обратной задачи в каждом из слоев будем использовать метод конечных разностей. На области изменения аргументов введем равномерную дискретную сетку
X, Уj, 2т, ік): X- = 'К
0 <- ; Уі = А 0 < і < N;
z т = тр, 0 < т < ;
їк = кг, 0 < к < М
где к - шаг сетки по горизонтальным координатам, р - шаг сетки по вертикальной координате, т - шаг сетки по времени. Параметры к и т связаны условием устойчивости Куранта с~05г < к . Параметры М, N задают количество точек сетки по горизонтальным координатам, М — количество слоев регистрации данных по вертикали.
В области, не содержащей источников, получаем разностную схему для дифференциального уравнения (1):
Vijmk
uijk+1 _ 2uijk + uijk_1
ui+1 jk _ 2uijk + ui_1 jk
uij+1k _ 2uijk + uij_1k h2
_ 0.
Здесь щк - значения и (г ^ )в точке (х, у) в момент времени ^, Су - значения с (г) в точке (хг, у). Расчет распространения звуковой волны (расчет «в прямом времени») выполняется по временным слоям в явной форме:
uijk+1 _ uijk(2 _
2х2 2х2
) +
(ui+1 jk + ui _1 jk)X 2 (ui.+1k + uj_1k)X 2
+------- ----ТГ1--------+ —1---------TT-----------иФ _1,
сЪ сМ2 j
(3)
начальные условия задаются в виде иг}-0 = = 0, граничные условия задаются в виде
(ui
+1jk ui_1jk
2h
для i = 1 или Nx-1,
) _ + (uijk+1 _ uijk _1)
с 052x
(uij+1k _ uij_1k) _ + (uijk+1 _ uijk_1)
2h
для j = 1 или Ny-1 .
с 052x
с
2
2
h
X
Рис. 1. Метод распараллеливания расчетов
Расчетная модель, описываемая формулой (3), аппроксимирует гиперболическое уравнение второго порядка (волновое уравнение) со вторым порядком точности. Эта модель учитывает волновые эффекты дифракции, рефракции, переотражения и т.д. Для решения задачи в каждом томографическом слое необходимо выполнить 0( №хМуМ1) операций. Таким образом,
объем вычислений растет не более чем линейно от числа точек сетки по времени N и числа точек сетки области вычислений К хЖу. Невязка на текущей итерации вычисляется по формуле
1 т
Ке*=2 2 2 (“«• - и,к )’^
к=0 (і,і)єХ
Здесь 5 - граница, - значения Щг^) в точке
(хг, у) на границе 5 в момент времени tk.
Зондирующий импульс задается формулой
И..0 = «п
^ 2п (Rіі-BegImp) WidtкImp
Л
- 0.5п
+1,
где Rij = (/ - /5ои, )2 + у - ЛБои,)2 ; параметр
WidthImp - длина волны излучения; /Беи,,,^ои, - координаты положения источников на границе 5, 5=1, ..., К5; Ж, — количество положений источника; Beglmp — расстояние в начальный момент времени от переднего фронта импульса до источника импульса.
Рассмотренный дифференциальный подход имеет ряд существенных преимуществ по сравнению с широко известным интегральным подходом, прежде всего, в объемах вычислений [14]. Использование, например, явных схем приводит к необходимости выполнить порядка 0(п4) операций типа сложения-умножения, где
п - размер области по одной координате. Эта оценка для числа операций на несколько порядков меньше, чем в интегральном подходе.
Исследование эффективности и масштабируемости программы
Возможность эффективного использования параллельных вычислений при решении обратной задачи определяется структурой алгоритма, возможностью выделения большого числа максимально независимых вычислительно блоков. Предложенная архитектура программы, позволила выделить большое количество таких блоков, на которых параллельно выполняются процессы общего потока задачи.
Во-первых, трёхмерная область делится на N независимо обрабатываемых слоев (М ~ 40) в соответствии с рис. 1.
Далее, в каждом слое расчёт для каждого источника (их количество К ~ 10-20) может быть выполнен в значительной степени независимо (результаты расчетов должны быть лишь просуммированы для вычисления градиента на каждой итерации градиентного спуска). Таким образом, в целом имеем набор М*К матриц, вычисление которых может быть выполнено независимо.
Для дальнейшего распараллеливания вычислений воспользуемся известным методом распараллеливания по пространству, который состоит в том, что общее поле вычислений размером Мх*Му точек матрицы разбивается на МуРаг*ШРаг=М частей-блоков размером МхЫ*МуЫ, вычисления в которых производятся различными вычислительными ядрами параллельно (М - число ядер) (рис. 1).
Один из основных принципов разбиения на блоки состоит в выравнивании загрузки
Таблица 1
Сравнение вариантов разбиения сетки, размер сетки - 502x502, _______________________число процессов 36 _______________________
Разбиение матрицы 1*36 2*18 3*12 6*6
Время, с 83 59 55 49
Таблица 2
Сравнение вариантов разбиения сетки, размер сетки - 1002x1002, _________________________число процессов 64 ___________________________
Разбиение матрицы 1*64 2*32 4*16 8*8
Время, с 692 386 254 220
процессоров. Для минимизации времени ожидания процессором данных с соседних процессоров расчёты всех блоков должны быть завершены синхронно, поэтому следует разбивать матрицу на блоки с равным объёмом вычислений. В нашем случае явной разностной схемы, одинаковой во всех точках, используем простейший вариант разбиения на блоки одинакового размера в соответствии с рис. 1.
Для минимизации межпроцессорных обменов следует минимизировать периметр блоков (в данном случае периметр равен 2*(МЫ + МЫ1) по отношению к их площади (в данном случае площадь равна МуЫ * МЫ). Как известно, минимум периметра у прямоугольника заданной площади достигается при МуЫ1 = МхЫ1. Это требование было проверено экспериментально на модельных задачах. Для этого были проведены эксперименты на сетках размером 502x502 и 1002x1002, результаты которых приведены в табл. 1 и 2. В верхней строке указаны значения МХ_РАЯ*Ш'_РАЯ для различных разбиений матрицы по процессам, в нижней - время расчета 10 итераций. При перестановке значений МХ_РАЯ и МУРАЯ результат не меняется. Из данных таблиц видно, что чем ближе значения МХ_РАЯ и МУ_РАЯ друг к другу, тем меньше время счета.
Для явной 7-точечной 4-связной в плоскости (х,у) разностной схемы при разбиении на прямоугольные параллельно обрабатываемые блоки возникает 4-связная сетка межпроцессорных обменов. На каждом процессоре (вычислительном ядре) для одного шага разностной схемы по времени выполняется ~МхЫ]1 арифметических операций и производится передача и приём ~МхЫ1 данных в 4 направлениях. Для обеспечения масштабируемости коммуникационная сеть кластера должна содержать, как минимум, такое же число связей, тогда все передачи данных могут выполняться параллельно.
Отметим, что поскольку в явной разностной схеме значения во всех точках сетки на следующем шаге по времени не зависят друг от друга и могут быть вычислены параллельно как SIMD-алгоритм, то задача эффективно выполняется на широком диапазоне архитектур про-
цессоров: общего назначения, векторных, массивно-параллельных, графических и т.п. Пример эффективного распараллеливания в задаче моделирования распространения волн от точечного источника для более сложной физической модели трехмерной упругой среды приведен в работе [18].
Рассмотрим вопросы эффективности вычислительной системы в задаче ультразвуковой томографии в зависимости от количества вычислительных ядер Nproc на одну матрицу (или в зависимости от размера блока Nxbl). Качественно опишем математической параметрической моделью процесс вычислений для нашей задачи.
Исходными параметрами, определяющими конфигурацию системы, являются:
• V - объём данных одного блока, V = Nxbl * * Nybl= Nx*Ny / NProc,
• Vc - объём кэша вычислительного узла;
• Vm - объём данных, не поместившихся в кэш, Vm = V - min(V,Vc);
• Bc - производительность вычислений в кэше (cache bandwidth);
• Bm - производительность вычислений в памяти (memory bandwidth) ;
• Lcomm - задержка распространения сообщений (comm latency);
• Bcomm - производительность коммуникационной сети (comm bandwidth).
Времена выполнения одного шага явной разностной схемы для одного процесса с точностью до коэффициентов равны:
• Tcache = min(V,Vc) / Bc — для операций в кэше;
• Tmem = Vm /Bm — для операций в памяти;
• Tcomm = 4 * Nxbl / Bcomm — для коммуникаций.
Общее время выполнения вычислений имеет вид
T= (Lcomm+Tcomm+Tcache+Tmem) x Nt, (4) где Nt - число шагов по времени разностной схемы.
Рассмотрим эффективность E функционирования каждого процесса в выбранной схеме
Таблица 3
Результаты тестирования на суперкомпьютере «Ломоносов»
Число МрЮе) процессов 1 2*2 3*3 4*4 5*5 6*6 7*7 8*8 9*9 10*10 11*11 12*12 13*13 14*14
Время (с) расчетов 15 итераций 1893 903 370 291 162 121 67 64 44 40 41 42 44 45
Рис. 2. График зависимости эффективности Е от Рис. 3. Зависимость производительности системы от Кргос на одной матрице: теоретическая зависимость числа процессоров на одну матрицу (штрих) (сплошная линия) и экспериментальная зависимость (штриховая линия)
распараллеливания. С точностью до коэффициента эффективность можно понимать как отношение времени То расчетов задачи на Ыргос процессах с пиковой производительностью к времени Т расчета на Ыргос процессах в выбранной схеме распараллеливания. Типично пиковой производительностью системы является производительность вычислений в кэш-памяти:
Е =То = К , (5)
Т Кргос • Т ^
где К - некоторый постоянный коэффициент, Т вычисляется по формуле (4) (зависит от Ыргос).
На рис. 2 приведено качественное поведение графика теоретической зависимости эффективности Е от Ыргос, построенное с использованием формул (4), (5) (сплошная линия), размер сетки матрицы Кх=Ыу = 1002. При построении графика в модели использованы типичные значения параметров для процессоров общего назначения.
Ниже в табл. 3 приведены результаты тестирования разработанной программы в рассматриваемой задаче ультразвуковой томографии, проведённые на суперкомпьютере «Ломоносов».
На основании данных, приведенных в табл. 3, на рис. 2 штриховой линией показан экспериментально полученный график эффективности в зависимости от числа процессов. Видно, что графики теоретической и экспериментальной кривых имеют похожие закономерности.
Первоначально при увеличении количества ядер на слой эффективность несколько падает,
поскольку при малом количестве ядер объем обменов по границам блоков невелик (для 1 ядра на матрицу обмены отсутствуют вообще). Далее эффективность начинает расти, что связано с увеличением доли вычислений в кэше. При дальнейшем увеличении числа процессов эффективность падает из-за увеличения ожидания процессом данных с соседних процессов.
Зависимость эффективности вычислений от количества процессов на слой зависит от параметров конкретной вычислительной системы. Тем не менее, существует некоторое пороговое значение количества процессов на матрицу (в нашей программе примерно 100 процессов), выше которого эффективность падает. Этот пик эффективности связан с расчетами в кэше. На современных вычислительных системах производительность арифметического устройства и кэш-памяти во много (10-100) раз выше производительности памяти, и соответствие размера данных размеру кэша оказывает сильное влияние на эффективность [19].
Знание этого порогового значения количества процессов в рассматриваемой задаче особенно важно в том случае, если количество доступных процессоров ограничено, поскольку правильное распределение процессоров по матрицам влияет на производительность вычислительной системы. Например, если доступно 200 вычислительных узлов, то использование их всех для вычисления на одной матрице или ис-
4а
Ф
5а
6а 6б
Рис. 4-6 (а, б). Слева - модельные сечения 3D объекта, справа - восстановленные сечения 3D объекта
пользование менее 50 узлов на матрицу было бы не эффективно.
Для вычислительных систем, в которых возможно параллельное выполнение обмена данными и вычислений, использование этой возможности должно приводить к существенному повышению эффективности расчетов. Однако проведенные эксперименты с неблокирующими коммуникационными операциями не привели к сокращению времени расчетов, что требует дополнительных исследований.
На рис. 3 штриховой линией показан график функции 1/Т, где Т(Щтс) - время расчетов Щтс процессоров на одной матрице. Этот график описывает также производительность, или ускорение вычислительной системы. Как видно, производительность растет с постоянной скоростью и выше вплоть до 100 процессоров на матрицу. Далее производительность ухудшается.
Численные расчеты
Рассматриваемая программа решения прямых и обратных задач томографической диагностики тестировалась в конфигурации, обес-
печивающей одновременную работу 20480 процессов на суперкомпьютере. Экспериментальные исследования проводились на компьютерно-синтезированном 3Э объекте с модельными неоднородностями. На рис. 4-6 приведены три модельных сечения (слева) и результаты реконструкции (справа) функции скорости распространения ультразвуковой волны у(х, у, 2) в исследуемом объекте как функции от х, у для трех фиксированных 2 = 21 (/ = 1, ..., 40). Расстояние между сечениями 21 (/ = 1, ..., 40) по оси 2 составляло 5 мм. Потемнение в каждой точке рисунка пропорционально у(х, у, 2). Минимальный размер неоднородности 3 мм. Вариация скорости у(х, у, 2) не превышала 20%. В ходе экспериментальных исследований программы решалась прямая задача распространения ультразвуковой волны в каждом слое. По полученным данным без внесения дополнительного шума решалась обратная задача восстановления функции у(х, у, 2) в каждом слое. Расчет прямой и обратной задач производился одновременно на 40 слоях. Параметры расчетной модели: длина волны излучения 5.0 мм; шаг регистрации сигналов по пространству 2.4 мм; размер облас-
ти ультразвукового зондирования по горизонтали 200x200 мм, по вертикали 200 мм.
Для решения обратной задачи использовался итерационный процесс с начального приближения v(x, y, z) = const - известной скорости вне области неоднородности. Количество итераций 700, время расчета около 4 часов. За 700 итераций функционал невязки (2) уменьшился в 10000 раз. Расчеты проводились на сетке 1002x1002 точки в 40 слоях для 8 источников излучения. Распараллеливания по пространству по координатам X и Y состояло в том, что общее поле вычислений размером 1002x1002 точки разбивалось на NyPar * NxPar одинаковых частей (NyPar = 8, NxPar = 8), вычисления в которых производились различными вычислительными ядрами. В результате получили распараллеливание задачи на 40 * 8 * 8 * 8 = 20480 процессов. В соответствии с предыдущим разделом, такое разбиение по процессам позволяет достичь высокой эффективности. Относительная эффективность распараллеливания на 20480 процессов составила около 60% (если сравнивать ускорение расчетов по отношению к расчетам на 1 процессоре).
Заключение
Продемонстрирована возможность эффективного использования разработанного программного обеспечения для большого числа процессоров. Программа обеспечивает высокий уровень масштабирования для нескольких десятков тысяч процессов.
Максимальная производительность каждого ядра (и производительность всей системы при ограниченном количестве ядер) достигается при определённом значении количества процессов на матрицу Nproc = Nopt, которое может быть определено для конкретной системы. При отклонении Nproc от оптимального значения как в большую, так и в меньшую сторону эффективность будет падать. Исходя из этого следует сначала выбирать Nproc, а затем, учитывая имеющиеся ресурсы системы, выбирать число параллельно рассчитываемых на системе матриц.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 12-07-00304 а. Статья подготовлена по материалам доклада автора на 14-я международная суперкомпьютерная конференция «Научный сервис в сети Интернет: поиск новых решений» (2012; http://agora.guru.ru/abrau2012)
Список литературы
1. Li C., Duric N., Littrnp P., Huang L. In vivo breast soundspeed imaging with ultrasound tomography // Ultrasound Med. Biol. 2009. V. 35. № 10. P. 1615-1628.
2. Jink R., PeterHk I., Ruiter N., et al. Sound-Speed Image Reconstruction in Sparse-Aperture 3-D Ultrasound
Transmission Tomography // IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control. 2012. V. 59. № 2. P. 254-264.
3. Wiskin J., Borup D.T., Johnson S.A., and Berggren M. Non-linear inverse scattering: High resolution quantitative breast tissue tomography // J. Acoust. Soc. Am. 2012. V. 131. № 5. P. 3802-3813.
4. Овчинников С.Л., Романов С.Ю. Организация параллельных вычислений при решении обратной задачи волновой диагностики // Вычислительные методы и программирование. 2008. Т. 9. № 1. C. 338-345.
5. Гончарский А.В., Овчинников С.Л., Романов С.Ю. Об одной задаче волновой диагностики // Вестник МГУ. Сер. 15. Вычисл. матем. и киберн. 2010. № 1. С. 7-13.
6. Бакушинский А.Б., Козлов А.И., Кокурин М.Ю. Об одной обратной задаче для трехмерного волнового уравнения // Журн. вычисл. матем. и матем. физики. 2003. Т. 43. № 8. С. 1201-1209.
7. Backushinsky A., Goncharsky A., Romanov S., Seatzu S. On the identification of velocity in seismics and in acoustic sounding // Pubblicazioni dell’istituto di analisa globale e applicazioni, Serie «Problemi non ben posti ed inversi». 1994. № 71.
8. Головина С.Г., Романов С.Ю., Степанов В.В. Об одной обратной задаче сейсмики // Вестник МГУ. Сер. 15. Вычисл. матем. и киберн. 1994. № 4. С. 16-21.
9. Bakushinsky A.B., Kokurin M.Yu. & Kozlov A.I. On stable iterative methods of gradient type for the inverse medium scattering problem // Inverse Problems in Science and Engineering. 2005. V. 13. № 3. P. 203-218.
10. Natterer F., Wubbeling F. A propagation-backpropagation method for uzltrasound tomography // Inverse Problems. 1995. V. 11. № 6. P. 1225-1232.
11. Beilina L., Klibanov M.V. Approximate global convergence and adaptivity for coefficient inverse problems. N.Y.: Springer, 2012.
12. Гончарский А.В., Романов С.Ю. Супер-компьютерные технологии в разработке методов решения обратных задач в УЗИ-томографии // Вычислительные методы и программирование. 2012. Т. 13. № 1. С. 235-238.
13. Воеводин Вад.В., Овчинников С.Л., Романов С.Ю. Разработка высокоэффективных масштабируемых программ в задаче ультразвуковой томографии // Вычислительные методы и программирование: новые вычислительные технологии. 2012. Т. 13. № 1. С. 307-315.
14. Гончарский А.В., Романов С.Ю. О двух подходах к решению коэффициентных обратных задач для волновых уравнений // Журн. вычисл. матем. и матем. физики. 2012. Т. 52. № 2. С. 1-7.
15. Кокурин М.Ю. О редукции нелинейной обратной задачи для гиперболического уравнения на плоскости к линейному интегральному уравнению // Вычислительные методы и программирование: новые вычислительные технологии. 2009. Т. 10. № 1. С. 300-305.
16. Ватульян А.О. О вариационной постановке обратных коэффициентных задач для упругих тел // Доклады Академии наук. 2008. Т. 422. № 2. С. 182184.
17. Ватульян А.О., Явруян О.В., Богачев И.В. Идентификация упругих характеристик неоднородного по
толщине слоя // Акустический журн. 2011. Т. 57. № 6. С. 723-730.
18. Глинский Б.М., Караваев Д.А., Ковалевский В.В., Мартынов В.Н. Численное моделирование и экспериментальные исследования грязевого вулкана «Г ора Карабетова» вибросейсмическими методами //
Вычислительные методы и программирование: новые вычислительные технологии. 2010. Т. 11. № 1. С. 95-104.
19. Foster I. Designing and building parallel programs: concepts and tools for parallel software engineering. Reading: Addison Wesley, 1995.
SCALABILITY OF THE PROGRAM FOR SOLVING THE INVERSE PROBLEM OF WAVE TOMOGRAPHY
S.Yu. Romanov
The aim of this work is to develop algorithms and scalable software for petaflop supercomputers to be used for solving nonlinear coefficient inverse problems of wave tomography with extremely large data volumes and super-high resolution. We propose a program structure that allows computations to be parallelized with a high degree of scalability. Some issues of computational efficiency are discussed and the results of our computations performed on 20480 computer cores of the «Lomonosov» supercomputer of the Supercomputing Center of Lomonosov Moscow State University are reported.
Keywords: coefficient inverse problems, wave equation, ultrasonic tomography, parallel computing, scalability, supercomputer.
