CC BY
7
Вычислительные технологии
Том 21, № 3, 2016
Решение задачи моделирования распространения возмущения в массивной породе на реконфигурируемой вычислительной системе
А. Е. Чистяков1, Ю.И. ДорончЕнко2, Д. А. Сорокин2'*
1НИИ многопроцессорных вычислительных систем ЮФУ, Таганрог, Россия 2Научно-исследовательский центр супер-ЭВМ и нейрокомпьютеров, Таганрог, Россия *Контактный e-mail: jotun@inbox.ru
Рассмотрен один из вариантов преодоления проблемы моделирования геофизических процессов — применение реконфигурируемых вычислительных систем (РВС), так как разработка новых месторождений в труднодоступных регионах требует быстрого и качественного анализа получаемых данных, что затрудняется необходимостью охвата широкого спектра геологических особенностей и учета их параметров в создаваемой геофизической модели. Данные вычислительные системы за счет адаптации архитектуры под структуру решаемой задачи позволяют получать существенно большую производительность вычислений, что особенно характерно для сильносвязанных задач, к числу которых относятся многие геолого-геофизические задачи.
Предложена новая реализация задачи моделирования распространения возмущений в массивной породе с различными неоднородностями, позволяющая эффективно масштабировать решение при увеличении доступного аппаратного ресурса РВС, использование которых позволяет значительно сократить время вычислений при сохранении точности результатов решения задачи. Последнее обеспечивает более быструю и качественную интерпретацию геофизических данных.
Ключевые слова: реконфигурируемая вычислительная система, моделирование волновых возмущений в неоднородной массивной породе, система линейных алгебраических уравнений, модифицированный попеременно-треугольный метод.
Введение
Большинство перспективных месторождений нефти и газа в Российской Федерации находится в труднодоступных регионах, что осложняет и значительно увеличивает во времени процесс их разработки. Одной из главных проблем, с которой при этом сталкиваются исследователи, является скорость интерпретации геофизических данных. При геолого-сейсмическом моделировании для труднодоступных регионов охватывается существенно более широкий спектр геологических особенностей. Поэтому скорость построения модели структурных построений границ раздела слоев в осадочной толще массивных пород и точность получаемых расчетов зачастую не удовлетворяют требованиям, предъявляемым интерпретаторами, поскольку перед ними стоит задача более качественного истолкования наблюдаемой картины. Для преодоления данной проблемы
© ИВТ СО РАН, 2016
необходимо совершенствовать методы и средства решения подобных задач в сейсморазведке, поскольку они возникают как на стадиях создания и детализации модели толщи массивных пород, так и на этапах верификации результатов физических наблюдений.
Одним из возможных вариантов повышения скорости решения задачи моделирования структурных построений границ раздела слоев при сохранении требуемой точности вычислений является применение реконфигурируемых высокопроизводительных вычислительных систем (РВС) на основе реконфигурируемых программируемых логических интегральных схем (ПЛИС). Главным преимуществом РВС перед традиционно применяемыми в данной области кластерными системами на основе процессоров общего назначения или графических ускорителей является возможность адаптации собственной архитектуры к структуре реализуемого алгоритма. За рубежом в последнее десятилетие активно ведутся исследования по созданию ПЛИС-ориентированных вычислительных систем, которые демонстрируют существенно большую эффективность при решении подобного класса задач.
Целью работы является сокращение времени вычислений и сохранение точности результатов решения задачи моделирования распространения звуковых возмущений в массивной породе, содержащей различные неоднородности, за счет применения ре-конфигурируемой вычислительной системы. Построена дискретная математическая модель распространения волновых геофизических процессов. Для повышения реальной точности решений использованы сетки, учитывающие заполненность расчетных ячеек. Дискретная модель создана на основе метода баланса. Предлагаемый численный метод имеет большую общность и пригоден для исследования процессов взаимодействия сейсмических волн с неоднородными включениями, поскольку позволяет наиболее корректно конструировать вычислительные алгоритмы на границах области интегрирования и раздела сред.
1. Постановка задачи
Для построения дискретной математической модели распространения волновых процессов в массивных породах требуется найти решение неоднородного волнового уравнения [1-3]
4 = ^2дс + /, (1)
удовлетворяющего начальным с(х, у, 0) = <р0(х, у) и граничным с(х, у, 0) = <р0(х, у) условиям при (х, у) Е 7. Здесь V — скорость распространения волновых колебаний; / — некоторая функция-источник волновых процессов; и — функции, определяющие соответственно с и с' в начальный момент времени £ = 0; 7 — функция, определяющая границы расчетной области.
Наиболее эффективный подход в решении подобного вида задач — это интегро-интерполяционный расчет в равномерной прямоугольной области [4, 5], при котором задача сводится к нахождению корней сеточных уравнений. Однако, помимо требований к скорости построения модели распространения волновых процессов, необходимо обеспечивать заданную точность вычислений, существенно влияющую на качество интерпретации. Поэтому для решения уравнений предлагается использовать один из наиболее актуальных методов, способных обеспечить "гладкость" решения задачи [6], а именно модифицированный попеременно-треугольный (МПТ) метод. Важной особенностью МПТ-метода является способность давать решение за наименьшее число итераций в случае плохо обусловленных задач, к которым относится и решаемая задача.
2. Суть МПТ-метода решения сеточных уравнений
В конечномерном гильбертовом пространстве Н рассматривается задача об отыскании решения операторного уравнения [6]
Ах = А : Н ^ Н, (2)
где А — линейный, самосопряженный (А = А*), положительно определенный оператор (А > 0). Для решения задачи (2) используется неявный итерационный процесс
В--— + Ахп = /, В : Н ^ Н, (3)
Гп+1
где В — некоторый обратимый оператор; п — номер итерации; тп+\ — итерационный параметр.
Оператор В рассчитывается по формуле
В = (Б + шп Я1)0-1(0 + ип Иа), (4)
где и > 0, у <Е Н.
Непосредственно алгоритм МПТ-метода (далее по тексту — алгоритм) для расчета сеточных уравнений описан в работе [7]. Рассмотрим вычислительную структуру данной задачи.
3. Вычислительная структура задачи решения сеточных уравнений МПТ-методом
Информационный граф задачи [8] состоит из n базовых подграфов П G = Уп, а каждый базовый подграф Пп можно разделить на шесть последовательно выполняемых подзадач (рис. 1).
Рассмотрим подзадачу Р1 формирования вектора невязки r = f — Ахп на основании коэффициентов сеточного уравнения A, B, Сп и F. Также здесь вычисляется равномерная норма невязки max.
Расчет элементов вектора невязки представляет собой линейный процесс
ТmO fmO ^тО^тО + ^1т0^т1 + ^2т0+ ^3т0+ ^4т0^т4 ,
в котором отсутствуют обратные зависимости. Это связано с тем, что rm0 = f (rm0±step) V step G {0,1, 2, 3,4,...}. Поэтому при структурной реализации подзадача Р1 пред-
Рис. 1. Базовый подграф задачи
ставляется в виде линейного вычислительного конвейера, структура которого показана на рис. 2. Здесь и далее через < гто > обозначен линейный поток параметров гто Е гп, где ш0 — адрес чтения, правила формирования которого описаны в [7]. Аналогично для
< /то >, где /то Е ^, < ато >, где ато Е А, < Ь1то > — < 64то >, где Ь1то - Ь4то Е В,
< Сто > — < Ст4 >, где Сто — Ст4 Е С.
Для расчета очередного значения гто Е гп выполняется 11 операций в формате РР32 (32-разрядной плавающей запятой стандарта 1ЕЕЕ754), что обусловлено требуемой точностью вычислений.
Из рис. 2 видно, что для конвейера подзадачи Р1 требуется 11 каналов доступа к данным (по одному каналу доступа к массивам А и четыре канала доступа к массиву В, три канала доступа к массиву С и два выходных канала < гто > и тах). В вычислительной структуре Р\ присутствуют блоки задержек 5 =12 и 5 =19, обеспечивающие выравнивание потоков данных на требуемое число тактов между вычислительными блоками конвейера.
Теперь перейдем к подзадаче Р2, которая формирует вспомогательный вектор уп алгоритма на основании коэффициентов А, В сеточного уравнения (И + шп)у = г, параметра шп и вектора невязки гп.
Расчет элементов вектора уп представляет собой рекурсивный процесс
У то = (Г то + (Ь2то Ут2 + 64то^т4)^)(«то (1 + ^/2))
-1
Это связано с тем, что уто = /(уто-1,уто-их). При структурной реализации построение линейного вычислительного конвейера невозможно, поэтому данная подзадача может быть реализована в виде вложенного конвейера [9, 10], структура которого показана на рис. 3.
Для реализации вложенного конвейера вокруг вычислительной структуры замыкается обратная связь, по которой циркулируют данные, пока не выполнятся все циклы алгоритма. При этом количество данных, циркулирующих по цепи обратной связи, равно глубине конвейера, данные подаются в вычислительную структуру порциями. После выполнения всех циклов формируется следующая порция данных для вычислительной структуры. Переключение каналов обратной связи осуществляется с помощью простейшего коммутатора.
Рис. 2. Структура вычислительного конвейера подзадачи
Особенность предлагаемой организации вычислений состоит в том, что данные во вложенный конвейер должны подаваться не порциями, равными глубине конвейера, а с изменяемой с каждой итерацией скважностью.
Латентность приведенного на рис. 3 вложенного конвейера определяется латентнос-тью блока ВР2 и составляет 92 такта. Поскольку скважность подачи данных на вход вычислительной структуры вложенного конвейера по мере его заполнения каждые 92 такта будет уменьшаться на 1, в начальный момент времени работы вычислительной структуры заполнение вложенного конвейера будет происходить 4278 такта. Поэтому входной поток <гт0 > необходимо буферизировать в памяти Ъи^ глубина которой должна быть, соответственно, не менее чем 4278 32-разрядных слов. После того как вложенный конвейер будет наполнен, данные на его вход можно подавать со скважностью, равной 1. Другими словами, на одной вычислительной структуре подзадачи Р2 за 92 такта будут выполняться 92 шага обработки данных. Структура блока ВР2 показана на рис. 4.
При расчете очередного значения ут0 Е уп выполняются 10 операций в формате ЕР32, а также требуются семь каналов-источников данных и один канал-приемник данных. Здесь стоит отметить, что подзадачи Р3, Р5 и Р6 по способу организации вычислений представляют собой абсолютно линейные вычислительные конвейеры, подобные Р\, а подзадача Р4 структурно полностью соответствует подзадаче Р2, поэтому детальный разбор их особенностей оставим за рамками настоящей статьи.
Дальнейший анализ показал, что организация вычислений в подзадаче Р2 по схеме "вложенный конвейер" позволяет выполнить объединение вычислительных структур
Рис. 3. Структура вложенного конвейера подзадачи
Рис. 4. Структура вложенного конвейера подзадачи
Рис. 5. Конвейеризированная вычислительная структура базового подграфа
подзадач Р1, Р2 и Р3 в единый вычислительный конвейер Р123. Также для подзадач Р4, Р5 и Р6 возможно объединение в Р456.
Учитывая вышеизложенное, получим конвейеризированную вычислительную структуру базового подграфа Пп, показанную на рис. 5.
Чтобы начать вычисления в конвейере Р456, необходимо закончить вычисления в Р123, поскольку правила адресации к массивам данных различаются, а именно для Р123 г = 1,NX — 1,j = 1,Ny — 1, а для Р456 i = Nx — 1,1, j = Ny — 1,1. Налицо непреодолимый разрыв в конвейерных вычислениях для одной задачи G, поэтому, чтобы не простаивало оборудование одного из конвейеров ( Р123 или Р456), на представленной вычислительной структуре решаются последовательно две задачи — G1(A1 ,В1 ,F1 ,С\,cj1 ) и G2(A2,В2,F2,С2п).
4. Результаты экспериментальных исследований
Оценка показала, что суммарное количество одновременно выполняемых операций в конвейеризированной вычислительной структуре Пп, обрабатывающих данные в формате FP32, равно 55, при этом требуется 12 каналов доступа к внешней памяти с суммарной пропускной способностью не менее 154 • 109 бит/с.
Производительность вычислительной структуры решения задачи можно оценить по формуле
Р(Пп) = Nop(nn) • vk,
где Nop (Пп) — количество блоков операций в вычислительной структуре Пп; ь> — тактовая частота работы вычислительной структуры; k — коэффициент поправки, учитывающий накладные расходы на обмен данными с динамической памятью, наполнение конвейера и др.
Экспериментальные исследования были проведены на действующей РВС "Тайгета", построенной на основе ПЛИС FPGA XC7VX485T фирмы Xilinx. Вычислительный модуль (ВМ) "Тайгета" выполнен в 19" форм-факторе высотой 2U и предназначен для создания высокопроизводительных многостоечных РВС. В состав ВМ "Тайгета" входят четыре платы вычислительного модуля (ПВМ), встроенная управляющая ЭВМ, а также подсистемы питания, управления, охлаждения и др.
Плата вычислительного модуля представляет собой 20-слойную печатную плату с двухсторонним монтажом элементов, на которой располагаются восемь ПЛИС типа XC7VX485T-1FFG1761, содержащих по 48.5 млн эквивалентных вентилей, микросхемы распределенной памяти SDRAM типа DDR2 общим объемом 4.5 Гбайт, интерфейсы LVDS и Ethernet и другие компоненты. На рис. 6 показаны ВМ и ПВМ "Тайгета".
Рис. 6. Вычислительный модуль (а) и ПВМ "Тайгета" (б)
Тактовая частота реализации составила у = 400 МГц, поправка на накладные расходы к = 0.85. Производительность на ВМ "Тайгета", при которой на одной ПЛИС ХС7УХ485Т была реализована одна вычислительная структура Пп, составила 32Р(Пп) = 598.4 Гфлопс.
Стоит отметить, что проведенный анализ показал возможность реализовать до 1000 устройств, выполняющих математические операции в формате ЕР32, на вычислительном ресурсе одной ПЛИС ХС7УХ485Т. Это более чем в девять раз превосходит требуемый ресурс для реализации одной вычислительной структуры Пп. Другими словами, имеется существенный потенциал для увеличения производительности за счет распараллеливания вычислений.
Поскольку одна вычислительная структура Пп требует 12 каналов памяти, при построении РВС с оптимальной архитектурой на базе ПЛИС ХС7УХ485Т необходимо организовать для ПЛИС максимально возможное количество каналов доступа к памяти. Данная ПЛИС имеет в своем составе 56 высокоскоростных приемопередатчиков СТХ с пропускной способностью от 8 Гбит/с и позволяет обеспечить потенциально возможную пропускную способность до 700 Гбит/с. Разумеется, это требует применения высокоскоростных микросхем памяти типа НМС [11], которые обеспечивают скорость обмена данными на уровне от 960 Гбит/с, что вполне удовлетворяет предъявляемым требованиям. Таким образом, на одной ПЛИС ХС7УХ485Т можно реализовать до четырех вычислительных структур Пп, что обеспечивает возможность создания РВС, решающей рассмотренную задачу с производительностью не менее 2240 Гфлопс.
Заключение
Предложена эффективная реализация на РВС задачи моделирования распространения возмущений в массивной породе, содержащей различные неоднородности, которая обеспечивает масштабируемость вычислительной структуры.
Полученные экспериментальные результаты позволяют сделать вывод о том, что при реализации рассмотренной задачи на РВС возможно добиться качественного скачка в скорости решения, сохраняя заданную точность вычислений. При этом выигрыш в скорости решения практически линейно возрастает при увеличении аппаратного ресурса РВС.
Применяя описанные технологии, можно построить РВС, ориентированную на решение целого класса геолого-геофизических задач, что, несомненно, будет способствовать более быстрой интерпретации наблюдаемых геофизических процессов при разработке перспективных месторождений полезных ископаемых.
Благодарности. Работа выполнена при финансовой поддержке проекта из бюджета Союзного государства в рамках реализации государственного контракта Министерства образования и науки Российской Федерации от 17 июня 2015 г. № 14.964.11.0001.
Список литературы / References
[1] Квасов И.Е., Панкратов С.А., Петров И.Б. Численное исследование динамических процессов в сплошной среде с трещиной, инициируемых приповерхностным возмущением, сеточно-характеристическим методом // Матем. моделирование. 2010. Т. 22, № 11. С.109-122.
Kvasov, I.E., Pankratov, S.A., Petrov, I.B. Computation of dynamic processes in continuous media with a crack initiated by the near-surface disturbance using grid-characteristic method // Matematicheskoe Modelirovanie. 2010. Vol. 22, No. 11. P. 109-122. (In Russ.)
[2] Владимиров В.С. Уравнения математической физики: Учебник для физ. и мех.-матем. спец. вузов. 4-е изд., испр. и доп. М.: Наука, 1981. 512 с.
Vladimirov, V.S. Equations of mathematical physics: Uchebnik dlya fiz. i mekh.-mathem. spets. vuzov. 4-e izd., ispr. i dop. Moscow: Nauka, 1981. 512 p. (In Russ.)
[3] Сухинов А.И., Зуев В.Н., Семенистый В.В. Уравнения математической физики. Таганрог: ТРТУ, 2009. 337 с.
Sukhinov, A.I., Zuev, V.N., Semenistyi, V.V. Equations of mathematical physics. Taganrog: TRTU, 2009. 337 p. (In Russ.)
[4] Самарский А.А. Теория разностных схем. М.: Наука, 1989. 432 с.
Samarskiy, A.A. Theory of difference schemes. Moscow: Nauka, 1989. 432 p. (In Russ.)
[5] Сухинов А.И., Чистяков А.Е., Шишеня А.В. Оценка погрешности решения уравнения диффузии на основе схем с весами // Матем. моделирование. 2013. Т. 25, № 11. С. 53-64.
Sukhinov, A.I., Chistyakov, A.E., Shishenya, A.V. Error estimation for solution of the diffusion equation based on the schemes with weights // Matematicheskoe Modelirovanie. 2013. Vol. 25, No. 11. P. 53-64. (In Russ.)
[6] Коновалов А.Н. Метод скорейшего спуска с адаптивным попеременно-треугольным пе-реобусловливателем // Дифференциальные уравнения. 2004. Т. 40, № 7. С. 953-963. Konovalov, A.N. Method of steepest decent with adaptive alternating triangular preconditioner // Differencial'nye Uravnenija. 2004. T. 40, No. 7. P. 953-963. (In Russ.)
[7] Сухинов А.И., Чистяков А.Е. Адаптивный модифицированный попеременно-треугольный итерационный метод для решения сеточных уравнений с несамосопряженным оператором // Матем. моделирование. 2012. Т. 24, № 1. С. 3-20.
Sukhinov, A.I., Chistyakov, A.E. Adaptive analog-SSOR iterative method for solving grid equations with nonselfadjoint operators // Matematicheskoe Modelirovanie. 2012. Vol. 24, No. 1. P. 3-20. (In Russ.)
[8] Каляев А.В., Левин И.И. Модульно-наращиваемые многопроцессорные системы со структурно-процедурной организацией вычислений. М.: Янус-К, 2003. 380 с.
Kalyaev, A.V., Levin, I.I. Modular-scalable multiprocessor system with structural-procedural organization of calculations. Moscow: Yanus-K, 2003. 380 p. (In Russ.)
[9] Каляев И.А., Левин И.И., Семерников Е.А., Шмойлов В.И. Реконфигурируемые мультиконвейерные вычислительные структуры. Ростов н/Д.: ЮНЦ РАН, 2008. 320 с. Kalyaev, I.A., Levin, I.I., Semernikov, E.A., Shmoilov, V.I. Reconfigurable multi-conveyer computer structures. Rostov-on-Don: Izd-vo YuNTs RAN, 2008. 320 p. (In Russ.)
[10] Доронченко Ю.И. Метод операционно-графового описания одновременных вычислений для многопроцессорных систем // Матер. Междунар. науч.-техн. конф. "Многопроцессорные вычисл. и управляющие системы — 2007". Таганрог: ТТИ ЮФУ, 2007. Т. 1. С. 11-17. Doronchenko, Yu.I. Method of operationally graph description of simultaneous computations for multi-processor systems // Intern. Sci.-Techn. Conf. "Multiprocessor Computer Systems — 2007". Taganrog: TTI YuFU, 2007. T. 1. P. 11-17. (In Russ.)
[11] http://www.hybridmemorycube.org/ (дата обращения 10.11.2015) http://www.hybridmemorycube.org/ (accessed 10.11.2015) (In Russ.)
Поступила в 'редакцию 24 февраля 2016 г.
Implementation of modelling for perturbation propagating in a massive rock on a reconfigurable computer system
CHISTYAKOV, ALEXANDR E.1, DORONCHENKO, YURIY I.2, SOROKIN, Dmitriy A.2'*
1 Scientific Research Institute of Multiprocessing Computing and Control Systems SFU, Taganrog, 347928, Russia
2Supercomputers and Neurocomputers Research Center, Taganrog, 347900, Russia * Corresponding author: Sorokin, Dmitriy A., e-mail: jotun@inbox.ru
The paper considers an alternative implementation for modelling of geophysical processes by means of reconfigurable computer systems (RCS), because development of new fields in hard-to-reach regions requires fast and qualitative analysis of input data. However, the developed geophysical model requires taking into account a wide range of geological peculiarities and their parameters, and becomes rather complex. Owing to adaptation of the architecture to the structure of the problem under consideration, RCSs provide a significantly higher performance, which is especially typical for tightly-coupled problems such as many geological and geophysical ones. We propose a new implementation for modelling of perturbation spreading in a massive rock with various heterogeneities, which allows effective scaling of the solution if the available RCS hardware resource is increased. Using RCSs we can significantly reduce the processing time, keeping calculation accuracy unchanged, and hence, faster and more qualitative interpretation of geophysical data can be done. In implementation for the problem of modelling, we applied the most effective approach based on integro-interpolation method within a uniform rectangular region. The approach provides the required rate of generation of the process for propagation of waves. The required accuracy of calculations, which considerably influences on the quality of interpretation, is provided by the modified alternating triangular method (ATM). An important feature of the ATM for providing smooth solution is capability to solve the problem during the less number
© ICT SB RAS, 2016
112
А. Е. HHCTSKOB, K. M. flopoH^eHKO, fl. A. CopoKUH
of iterations when the geological-geophysical computational problem is ill-conditioned. Experimental research was performed on the RCS Taygeta designed on the base of up-to-date Xilinx FPGAs. With the help of the RCS high processing rate and keeping the prescribed accuracy of calculations are possible. Most importantly, we have proved that application of such technologies for implementation to such problems is really promising.
Keywords: reconfigurable computer system, modeling of ripple disturbances in an inhomogeneous massive rock, system of linear algebraic equations, modified alternating triangular method.
Acknowledgements. The paper was financially supported by the project from the budget of the Union State under the contract by Russian Ministry of Education and Science from June 17, 2015, No. 14.964.11.0001.
Received 24 February 2016
