CC BY
126ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
ПРИМЕНЕНИЕ ТЕХНОЛОГИЙ ПАРАЛЛЕЛЬНОГО ПРОГРАММИРОВАНИЯ ДЛЯ МОДЕЛИРОВАНИЯ СЕЙСМИЧЕСКИХ ВОЛН СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИМ МЕТОДОМ
1Иванов А.М., ^Хохлов Н.И., 1,2Петров И.Б.
Московский физико-технический институт, http://mipt.ru г. Долгопрудный, 141701 Московская область, Российская Федерация
2Научно-исследовательский институт системных исследований, Российская академия наук, https://www.niisi.ru 117218 Москва, Российская Федерация
Поступила в редакцию 25.11.2016
В работе рассматривается решение гиперболической системы уравнений сеточно-характеристическим методом применительно к задаче моделирования распространения сейсмических волн в твердых деформируемых телах. На примере этой задачи демонстрируется применение технологий параллельного программирования для систем с общей памятью с использованием технологий OpenMP и POSIX Threads. Для достижения максимальной производительности использовались векторные инструкции SSE и AVX на центральных процессорах. Также было произведено распараллеливание для графических процессоров при помощи технологий CUDA и OpenCL. На графических процессорах было рассмотрено распараллеливание на нескольких GPU, в том числе и с использованием технологии NVIDIA GPUDirect, позволяющей графическим процессорам осуществлять обмены данными напрямую без участия центрального процессора. Показаны достигнутое ускорение распараллеливания и процент от теоретически возможной производительности. В работе описываются произведенные оптимизации, которые позволили достигнуть полученных результатов ускорения.
Ключевые слова: математическое моделирование, параллельное программирование, сеточно-характеристический метод, общая память, сейсмика, векторизация, GPU, CUDA, OpenCL
УДК 519.633
Содержание
1. Введение (185)
2. математическая модель (186)
3. Численный метод (186)
4. Постановка задачи (187)
5. условия тестирования распараллеливания (187)
6. Процесс оптимизации (188)
7. Параллельная версия (189)
8. оптимизация алгоритма на графических процессорах (189)
8.1. Перенос реализации c CPU на GPU -cudal (190)
8.2.Структурамассивовипоследовательное обращение к памяти - cuda2 (191)
8.3. Параметры вызова kernel'a - cuda3 (191)
8.4. Размеры блоков - cuda4 (192)
8.5. Использование opencl1 (192)
9. Несколько графических процессоров (192)
10. Заключение (193) Литература (194)
1. ВВЕДЕНИЕ
Задача распространения волн в линейно-упругой среде представляет интерес в сейсмологии и в разведочной геофизике. С ростом производительности современных вычислительных систем и возможности параллельного выполнения программного кода, появилась возможность проводить компьютерное моделирование сейсмических волн с большой точностью. Возник вопрос, как эффективно использовать доступные вычислительные ресурсы.
Одно из решений этого вопроса — использование технологий параллельного программирования для систем с общей памятью. В статье [1] рассматривается задача распараллеливания моделирования эффектов, происходящих в грунте в ходе землетрясения, с использованием технологии
OpenMP. Применительно к другим задачам математического моделирования, с помощью технологии OpenMP в работе [2] распараллеливается конечноэлементный метод на неструктурированных сетках.
Другой подход — распараллеливание на современных графических процессорах. Этому посвящено большое число работ. Например, в статье [3] рассматривается распараллеливание конечно-разностного метода на GPU с использованием технологии CUDA [4]. Есть работы [5, 6] в которых решаются задачи обратной миграции на кластере из графических процессоров.
И, наконец, третий подход распараллеливание в системах с распределенной памятью. В таких системах также могут использоваться как центральные процессоры, так и графические. В [7] распараллелен конечно-разностный метод моделирования вязко-упроугой среды с использованием технологии MPI [8], а в работе [9] — метод спектральных элементов для моделирования распространения волн в астероиде. Кроме того, в работе [10] сравнивается эффективность распараллеливания на GPU с использоавнием CPU кластеров. Исследование [11] демонстрирует использование гибридной системы OpenMP + MPI для моделирования сейсмики конечно-разностным методом.
Существуют также и другие подходы, подробное рассмотрение которых остается за рамками данной статьи. Например, использование FPGA [12] и ASIC [13].
В данной работе рассматривается решение гиперболической систем уравнений, которая описывает поведение волн в твердых упругих телах. Для нахождения решения этой системы применяется сеточно характеристический метод [14—16]. Вычисления производятся по явной схеме, поскольку в таком случае программная реализация хорошо поддается распараллеливанию. Решение задачи распараллелено в системах с общей памятью с помощью технологий OpenMP, POSIX Threads [17] и на графических процессорах с использованием CUDA, OpenCL [18].
2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Для описания поведения среды использовалась модель идеального изотропного линейно-упругого материала. Приведённая ниже система дифференциальных уравнений в частных производных описывает состояние элементарного объёма упругого материала в приближении малых деформаций для двумерного случая:
Р
dv да дах
ду,, да да.
dt дх
ду
Р
дt дх
ду
да ду дуу
dt дх ду
дауу ду дуу
= + (Л + 2ß)—-,
dt дх ду
да
(
ху
&
= ß
дух дУу Л
— + — , дх ду )
где р — плотность среды, X, ^ — параметры Ламе, V и V — горизонтальная и вертикальные составляющие скорости частиц среды, охх, о , о — компоненты тензора напряжения.
Данную систему можно переписать в матричной форме:
ди
ди„
_JL + A _
dt pq дх
+ B
ди„
= 0,
pq ду ' (1)
где u — вектор из 5 независимых переменных u
— (о , о , о , v , v )T. Явный вид матриц A и B
^ xx^ yy xy 3 x 7 W L ÜQ
pq
PQ
представлен в [19]. Здесь и далее подразумевается суммирование по повторяющимся индексам. Собственные значения матриц А и В таковы:
г pq pq
S--c ,
1 p'
S2 = - С
—
S —c.
4 s-
s_ — c , где c
5 p p
и с — скорости распространения продольных и поперечных волн в среде.
3. ЧИСЛЕННЫЙ МЕТОД
Применяя покоординатное расщепление, можно свести задачу построения разностной схемы для системы уравнений (1), к задаче построения разностной схемы для системы вида:
ди
+ A _
dt pq дх
= 0,
J2)
Для гиперболической системы уравнений (2) матрицу A можно представить в виде A = RAR-1, где Л — диагональная матрица, элементы которой — собственные значения матрицы A, а R — матрица, состоящая из правых собственных
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
векторов матрицы А. Введём новые переменные: = К-1и (так называемые инварианты Римана). Тогда система уравнений (2) сведётся к системе из 5 независимых скалярных уравнений переноса.
Приведем схему третьего порядка точности для численного решения одномерного линейного уравнения переноса и + аи= 0, а > 0, о = ат/Ь, т — шаг по времени, Ь — шаг по координате:
и„
= < + а(Л о + А 2)/2 + а2(А 0-А 2)/2 +
(3)
(Л>-2Ло +Л2),
6
Л = и" , - и",
0 т—1 т'
Лп п
, = и . — и ,,
1 т— 2 т—1 >
Л = и" — и" + .
2 т т+1
Схема (3) устойчива для чисел Куранта, не превышающих единицу. Используется сеточно-характеристический критерий монотонности, опирающийся на характеристическое свойство точного решения:
шт«—1, ит) < ит+1 < тах(ит—1, ипт).
В местах выполнения данного критерия порядок схемы падает до второго.
После того как значения инвариантов Римана на следующем шаге по времени найдены,
востанавливается решение: ип
= К"№п
4. ПОСТАНОВКА ЗАДАЧИ
Тестовая модель изображена на рис. 1. Размеры области даны в километрах. На нижней и боковых границах устанавливалось неотражающее граничное условие, на верхней — свободная граница. В качестве источника возмущения использовалась вертикальная сила, приложенная к площадке с 925.7 м по 974.1 м на дневной поверхности, амплитуда которой описывалась импульсом Рикера частоты 40 Гц. Результаты расчета представлены на рис. 2.
5. УСЛОВИЯ ТЕСТИРОВАНИЯ РАСПАРАЛЛЕЛИВАНИЯ
Все тесты на центральном процессоре производились на одной и той же тестовой задаче. При этом использовалась двумерная сетка размерами 4000 х 4000 узлов. Выполнялось 100 шагов по времени. Отдельно были измерены результаты для вычислений с одинарной и с двойной точностью. На графиках в данной работе представлены резульаты для двойной точности расчетов. Каждый узел содержал по 5 переменных, поэтому в случае вычислений с одинарной точностью вся сетка занимала в памяти 305 Мбайт, а в случае двойной точности - 610 Мбайт.
Использовались следующие флаги компиляции: -Йпэ-^ее-уейэпге, чтобы запретить векторизацию там, где мы ее не хотим использовать, -бэрептр и -pthread для
Рис. 1. Геологическая модель антиклинальной ловушки [20]
— волновая картина в момент / = 038 с.
Таблица 1
Центральные процессоры, на которых производилось распараллеливание
Название (обозначение на графиках) Частота, ГГц Кеш, Кбайт Процессоры и ядра Компиляторы Производительность, Гфлопс
AMD Opteron 6272 (a64) 2.1 2048 8 процессоров по 8 ядер gcc4.4.6 1075.2
Intel Xeon E5-2697 (i24) 2.7 30720 2 процессора по 12 ядер gcc4.4.7, icc15.0.0 518.4
AMD Opteron 8431 (a48) 2.4 512 8 процессора по 6 ядер gcc4.4.7 921.6
задействования той или иной техологии, -O2 для использования стандартных оптимизаций, -msse, -mavx, чтобы явно указать какой набор инструкций мы хотим использовать при автовекторизации.
Измерения были произведены на процессорах, указанных в Табл. 1. В дальнейшем обозначения используемых компиляторов на графиках — gcc и icc. Обозначения для используемого набора векторных инструкций процессора — sse или avx.
Для тестов на GPU рассматривалась задача с другими параметрами. Это была двумерная задача с количеством узлов 4096^4096. Осуществлялось 6500 шагов по времени. В каждом узле сетки хранилось по 5 переменных с плавающей точкой. Все вычисления производились как с одинарной (SP), так и с двойной (DP) точностью. Размер сетки в памяти — 320 Мбайт для вычислений с одинарной точностью и 640 Мбайт для вычислений с двойной точностью. В Табл. 2 представлены характеристики графических процессоров, на которых производились тестовые расчеты.
Таблица 2
Характеристики используемых графических карт
Графический процессор Потоковые процессоры Тактовая частота, МГц ГФлопс (SP) SP:DP ГФлопс (DP)
GeForce GT 640 384 900 691 24 29
GeForce GTX 480 480 1401 1345 168
GeForce GTX 680 1536 1006 3090 24 129
GeForce GTX 760 1152 980 2258 24 94
GeForce GTX 780 2304 863 3977 24 166
GeForce GTX 780 Ti 2880 876 5046 24 210
GeForce GTX 980 2048 1126 4612 32 144
Tesla M2070 448 1150 1030 2 515
Tesla K40m 2880 745 4291 3 1430
Tesla K80 2496 562 2806 1.5 1870
Radeon HD 7950 1792 800 2867 4 717
Radeon R9 290 2560 947 4849 8 606
6. ПРОЦЕСС ОПТИМИЗАЦИИ
В первую очередь была написана последовательная версия алгоритма. Пересчет узлов производился в два шага: по оси X (обход сетки в начале по столбцам, потом по строкам) и по оси Y (обход сетки в начале по строкам, потом по столбцам). Было измеренно количество арифметических операций в коде программы и из этого было полученно теоретическое значение флопс — 262, требуемое для пересчета одного узла. На графиках эта версия обозначается simple.
Далее была предпринята попытка оптимизировать шаг по Y. На этом шаге был изменен порядок обращения к узлам так, что он совпал с обращением в шаге по X. Таким образом, обращение к памяти стало последовательным, что привело к более эффективному использованию кеш-памяти. Кроме того, путем замены арифметических операций на эквивалентные, было уменьшено требуемое количество флопс для пересчета узла до 190. Обозначение этой версии на графиках - cf.
После этого были векторизованы циклы, в которых обрабатывались узлы сетки. Это было сделано с помощью директив (#pragma omp simd) технологии OpenMP, появившихся в стандарте 4.0. Явно было указано какой набор инструкций должен использоваться: SSE или AVX. Измерения производились только на процессоре Intel Xeon E5-2697.
Исходя из частоты ядра каждого конкретного процессора была определена его пиковая производительность. Таким образом был вычислен процент от пиковой производительности последовательной версии для каждого теста (рис. 3).
a48, gcc, cf a48, gcc, simple a64, gcc, cf a64, gcc, simple i24, gcc, cf i24, icc, cf i24, gcc, simple i24, icc, simple i24, icc, cf, sse i24, icc, cf, avx
4 6 8 10 12 14 16
Рис. 3. Процент от пиковой производительности.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
ш s х ш
О ^
и >
20
10
9 -
•
/
•
A ■ -d —
й л i
/ / ri ■
J » r- ■
Л Л Ш- m
✓g ........... ■
/ -
и- — 1- ■ ч
1 г
■
20
30
40
л
ь
0
1
ш S
£ ш
m
05
Число потоков
Число потоков
-•- simple,gcc,aff -■- simple,gcc cf,gcc,aff * cf,gcc
-•- spl-gr,gcc . • spl-gr.gcc.aff -Ш- pthread,gcc
Рис. 4.
7. ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ
Алгоритм был распараллелен с использованием технологий OpenMP и POSIX Threads. При этом отдельно была распараллелена каждая из версий, чтобы была возможность пронаблюдать последствия оптимизации последовательного кода на эффективность его распараллеливания.
Была сделана оптимизация уже параллельного кода. В ней каждый поток выделяет память, необходимую только для пересчета узлов, отданных в распоряжение именно этого потока. На графиках эта версия обозначается spl-gr. Теоретически это должно привести к росту производительности в NUMA системах, так как память, выделенная конкретным потоком, согласно стандарту OpenMP выделяется в локальной памяти потока, обращение к которой занимает наименьшее количество времени.
на AMD Opteron 8431.
Отдельно были произведены измерения, в которых потокам было запрещено перемещаться между ядрами процессора, т.е. была осуществлена "привязка" потоков к определенным ядрам процессора. Это сделано для уменьшения временных затрат. Без этого оптимизация, связанная с разделением сетки на отдельные части для различных процессов перестает иметь смысл. Результаты тестов представлены на рис. 4-7.
Далее OpenMP версия кода с разделением сетки между потоками была изменена таким образом, чтобы использовать POSIX Threads в качестве технологии распараллеливания.
8. ОПТИМИЗАЦИЯ АЛГОРИТМА НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
За основу для реализации алгоритма на графических процессорах была взята оптимизированная для
выполнения
на
1
05
К
0 10
20
..........*
3
4
5
60
Число потоков
Число потоков
-•- simple,gcc,aff -■- simple,gcc cf,gcc,aff -^cf,gcc
spl-gr,gcc _spl-gr,gcc,aff pthread,gcc
Рис. 5. Результаты распараллеливания на AMD Opteron 6272.
1
0
0
0
0
1
0
ш s х ш œ о
U >
15
10
II!
Ц^Ш-Щ
il
tí
■û ь
0
1
со s
Ë ш
m
05
10 15 20
Число потоков
25
r
10 15 20
Число потоков
25
simple,gcc -■- cf,gcc -e-simple,icc cf,icc —♦— spl-gr,icc
cf,icc,sse -■- spl-gr,gcc • cf,îcc,avx_ * spl-gr,îcc,avx_ ♦ spl-gr,icc,sse
Рис. 6. Результаты распараллеливания на Intel Xeon E5-2697 (без привязки).
центральном процессоре версия данной программы. Оптимизациям была подвергнута наиболее затратная в вычислительном плане часть алгоритма. Так как использовалось расщепление по пространственной координате, для пересчета всей сетки требовалось два шага: по оси X и по оси Y. При этом было вычислено количество арифметических операций с плавающей точкой, требуемое для пересчета одного узла сетки за два шага — 190 Флопс. Таким образом, зная количество узлов сетки, количество шагов по времени, можно было определить теоретически потребляемое алгоритмом количество ГФлопс. Далее, зная количество потоковых процессоров в графическом процессоре, их тактовую частоту и число FMA (fused multiply-add) процессоров в одном процессоре, была вычислена пиковая производительность для каждого GPU. Реальные тесты алгоритма на графических процессорах показали меньшие значения
производительности. На рис. 8 приведен результат тестирования. Процент от пиковой производительности алгоритмов измерялся как соотношение между этими двумя значениями — теоретически требуемым количеством Флопс для пересчета сетки и реально потребленным количеством Флопс.
8.1. Перенос реализации c CPU на GPU -cudal
В первоначальном варианте алгоритм был перенесен на выполнение на графических процессорах с использованием технологии CUDA, но не оптимизирован для выполнения на GPU. В ней на графическом процессоре выделялось в 2 раза больше памяти, чем требовалось для хранения расчетной сетки. Это стандартная практика при работе с технологиями написания алгоритмов для графических процессоров. В результате происходит
05
Число потоков
Число потоков
-•-s¡mple,gcc • cf,¡cc,sse
cf,gcc -«-s¡mple,¡cc cf,¡cc —♦— spl-gr,icc
spl-gr,gcc • cf,¡cc,avx_ * spl-gr,¡cc,avx_ ♦ spl-gr,¡cc,sse
Рис. 7. Результаты распараллеливания на Intel Xeon E5-2697 (с привязкой).
1
5
0
0
0
0
5
5
1
0
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Radeon R9 290 Radeon HD 7950 Tesla K80 Tesla K40m Tesla M2070 GeForce GTX 980 GeForce GTX 780 Ti GeForce GTX 780 GeForce GTX 760 GeForce GTX 680 GeForce GTX 480 GeForce GT 640
Рис. 8.
cudal cuda2 cuda3 cuda4
opencll
10 20 30
' от пиковой m
синхронизация только между вызовами функций, исполняющихся на графическом процессоре (CUDA kernels). Это было сделано по причине того, что глобальная синхронизация всего графического процессора создает большие временные задержки, поэтому эти задержки были вставлены между вызовами kernel^. Таким образом, удалось уменьшить количество глобальных синхронизаций до двух раз на один шаг по времени. Архитектура потоковых процессоров CUDA такова, что в одном CUDA блоке потоки, исполняющие один и тот же код над разными областями памяти, выполняются одновременно, если нет ветвлений в коде программы. Поэтому ситуации, когда все потоки ждут завершения выполнения одного, случаются, только если этот один поток выполняет какие-то отличные от остальных операции.
Все операции по выделению памяти на GPU и вызовам функций, исполняющихся на графическом процессоре, производятся хостом — CPU. Операции копирования сетки из памяти хоста и обратно требуют много времени, поэтому сетка копируется один раз из памяти хоста в память графического процессора, перед началом основных вычислений, и один раз в конце из памяти GPU в память хоста. Количество шагов по времени и размер расчетной сетки были таковы, что время, требуемое для расчетов, намного превосходило время, требуемое на это копирование.
Данные в сетке можно хранить двумя способами: в виде массива структур и в виде структуры массивов. Под структурой в данном случае подразумевается вектор u, состоящий из 5 компонент. В первоначальной версии, в расчетной сетке данные были организованы в виде массива структур, то есть данные определенного узла лежат в памяти последовательно.
Для задания граничных условий задачи требовалось отличное от остальных потоков поведение на границе расчетной сетки. Для обработки границ сетки, в коде kernel^ помещен участок кода с пятью условными блоками if-else.
Для этой версии было перебрано несколько вариантов размера блока. Был найден оптимальный размер 16x16, при котором время работы данной версии было минимально.
8.2. Структура массивов и последовательное обращение к памяти — cuda2
Следующая оптимизация состояла в использовании общей памяти (shared memory) блока CUDA. Это делается по причине того, что задержка (latency) при взаимодействии с этой памятью меньше, чем при взаимодействии с глобальной памятью. Оптимизация заключается в том, что чтение из глобальной памяти происходит только один раз на каждом шаге по X и по Y в начале вызова функции, выполняющейся на GPU, при этом данные копируются в общую память. Сразу же после этого осуществляется синхронизация всех потоков в блоке, после чего все вычисления производятся над данными в общей памяти блока. В конце kernel^ результат записывается на вторую копию расчетной сетки в глобальной памяти.
Другая оптимизация заключалась в выборе другого способа хранения расчетной сетки в памяти графического процессора — структуры массивов. После этого обращения к глобальной памяти стали последовательными (coalesced), что привело к повышению производительности.
Указанные выше оптимизации позволили уменьшить количество блоков if-else при обработке границ сетки до двух.
8.3. Параметры вызова kernel'a - cuda3
До этого дополнительные данные, такие как размер сетки, материал, значения переменных, полученных в результате промежуточных вычислений, справедливых на протяжении всех шагов по времени, передавались в качестве указателя на структуру в глобальной памяти графического процессора. Все потоки постоянно обращались к одному и тому же участку памяти. В версии алгоритма CUDA3 такие данные вычислялись на CPU один раз перед началом
192 -
1J* ИВАНОВ а.М., ХОХЛОВ НИ., ПЕТРОВ И.Б. ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
выполнения кода на GPU и передавались через параметры вызова кегпеГа. Ожидалось, что у каждого потока будет создана локальная копия этих значений.
8.4. Размеры блоков - cuda4
Важно было подобрать размеры блоков так, чтобы графический процессор был постоянно загружен, то есть чтобы не было частично не занятых потоковых процессоров. Также важно было выбрать размеры блоков, кратные размерам warp^, так как внутри warp'a осуществляется последовательный доступ к памяти. Другая причина — потоки одного warp'a способны выполнять код одновременно над разными участками памяти, если нет ветвлений.
В версии CUDA4 подбираются размеры блоков таким образом, чтобы удовлетворить указанным выше требованием и минимизировать количество узлов, для которых требуются обмены памятью между блоками. В шаге по X для пересчета каждого узла расчетной сетки требуются значения в двух смежных узлах по оси X, а в шаге по оси Y — требуется два смежных узла по оси Y. Поэтому следовало подобрать размеры блоков так, чтобы количество узлов, для которых требуются значения в смежных блоках, было как можно меньшим. Это необходимо было сделать, так как для узлов смежных блоков требуется выделение дополнительной памяти в области общей памяти блока.
Если, например, взять блок размера M х N, то в шаге по X для хранения узлов смежных блоков потребуется 4N дополнительных узлов, а в шаге по Y — 4M . Поэтому в шаге по X размер был выбран равным 256 х 1, в результате блок требовал всего 4 дополнительных узла в памяти. В шаге по Y размер блока был принят равным 16X16, чтобы достичь компромиса между количеством дополнительной памяти (64 узла сетки) и требованиями к последовательному доступу к памяти.
8.5. Использование opencll
Следующим шагом было написание OpenCL реализации данного алгоритма. За основу для этого была взята оптимизированная версия алгоритма CUDA4.
Результаты тестирования ускорения для всех реализаций представленны на рис. 9.
Radeon R9 290 Radeon HD 7950 Tesla K80 Tesla K40m Tesla M2070 GeForce GTX 980 GeForce GTX 780 Ti GeForce GTX 780 GeForce GTX 760 GeForce GTX 680 GeForce GTX 480 GeForce GT 640
0
L cudal ^cuda2 ^cuda3 ^cuda4
•m opencll
Рис. 9. Ускорение на GPU по сравнению с CPU. Максимальное полученное ускорение по сравнению с CPU на одном графическом процессоре — в 55 раз на GeForce GTX 780 Ti в вычислениях с одинарной точностью и в 44 раза на Tesla K80 в вычислениях с двойной точностью. Следует отметить хорошие результаты для устройств от AMD, несмотря на то, что использовались дешевые десктопные карты, они показали хороший результат на реализациях и одинарной и двойной точностью.
9. НЕСКОЛЬКО ГРАФИЧЕСКИХ ПРОЦЕССОРОВ
Далее алгоритм был распараллелен для работы на нескольких графических процессорах. Для выполнения на нескольких GPU использовался наиболее оптимизированный вариант. Расчетная сетка при этом была разделена на несколько равных по размеру прямоугольных участков. Деление производилось по оси Y, вследствие выбранного размера блока по оси X.
Тесты с использованием нескольких GPU были сделаны только для одинаковых графических процессоров. Это сделано из-за того, что при использовании GPU с разной производительностью более быстрые процессоры окажутся причиной простоя более медленных графических процессоров, и эффект от их совместного использования может оказаться менее заметным. Синхронизация расчетной сетки между графическими процессорами производилась путем
обмена через память хоста (CPU). Причем синхронизация производилась только один раз на каждом временном шаге перед шагом по оси Y. Если размер сетки равен M x N, а число процессоров — D, то число требуемых для синхронизации узлов сетки — 4M(D — 1).
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Tesla K40m GeForce GTX 780 Ti Tesla M2070 GeForce GTX 680 Tesla K80 GeForce GTX 980
11 GPU □ 2 GPU L3 GPU 4 GPU 15 GPU
6 GPU
7 GPU
8 GPU
0
Рис. 10. Ускорение при использовании GPUDirect.
Также была реализована версия на основе GPUDirect. Основное преимущество технологии GPUDirect при решении поставленной задачи заключается в возможности пересылки данных, располагающихся в памяти графических процессоров, напрямую без участия хоста через шину PCI Express, т.е. нет необходимости копировать данные вначале с графического процессора на хост, а потом с хоста на другой графический процессор. Результаты тестирования приведены на рис. 10.
При использовании технологии CUDA результат отличается незначительно, при использовании в качестве буфера обмена памяти хоста. Для технологии OpenCL результат несколько хуже.
10. ЗАКЛЮЧЕНИЕ
В работе показано использование возможностей центральных процессоров и графических процессоров для решения задач сейсмики сеточно-характеристическим методом.
Наибольшая производительность от пиковой для процессора была получена в последовательной версии 22% для вычислений с одинарной точностью и 17% для вычислений с двойной точностью с оптимизацией доступа к кеш памяти и использованием векторных инструкций avx.
Максимальные полученные ускорения: на AMD Opteron 6272 - в 37 раз на 64 ядрах, на AMD Opteron 8431 - в 25 раз на 48 ядрах, и на Intel Xeon E5-2697 - в 17 раз на 24 ядрах.
Таким образом удалось продемонстрировать, что в системах с общей памятью основным ограничением для ускорения с ростом числа потоков является максимальная скорость взаимодействия с памятью. Разделение памяти между потоками позволяет не обращаться ядрам процессора к нелокальной памяти,
а значит доступ к памяти осуществляется быстрее. Привязка потоков к определенным ядрам уменьшает количество обращений в память, используемых для поддержания когерентности кешей ядер, и также ведет к увеличению производительности. Следует заметить, что реализация параллельной версии с использованием POSIX Threads оказывается сложнее, чем при использовании OpenMP, но при этом не дает значительного прироста производительности или эффективности распараллеливания.
В дополнение к этому, в работе были описаны методы, которые позволили достигнуть наибольшей производительности указанного алгоритма при вычислениях на GPU. Рассмотрен вопрос эффективного использования памяти графического процессора как в случае использования одного GPU, так и в случае использования нескольких графических процессоров. Было определено влияние различных оптимизаций на производительность алгоритма. Максимальное полученное ускорение по сравнению с CPU на одном графическом процессоре — в 55 раз на GeForce GTX 780 Ti в вычислениях с одинарной точностью и в 44 раза на Tesla K80 в вычислениях с двойной точностью. Для одинарной точности удалось достичь производительности в 460 ГФлопс, а для двойной точности была получена максимальная производительность в 138 ГФлопс. Максимальное достигнутое ускорение от количества графических процессоров — в 7.1 раз на 8 графических процессорах для двойной точности. Технология GPUDirect повысила ускорение до 10% от того, что было достигнуто без нее при вычислениях с одинарной точностью.
На основе полученных результатов можно сделать выводы о возможности применения графических процессоров для такого рода задач. Полученные результаты будут похожи для других задач, использующих похожие численные методы (метод конечных объемов, конечно-разностные методы). Также стоит отметить о хороших результатах для процессоров от производителя AMD и технологии OpenCL, в то время как большинство работ используют технологию CUDA для графических процессоров от NVidia.
2
4
6
Исследование выполнено при финансовой поддержке
РФФИ в рамках научного проекта № 15-07-01931 A.
ЛИТЕРАТУРА
1. Caserta A, Ruggiero V, Lanucara P. Numerical modelling of dynamical interaction between seismic radiation and near-surface geological structures: a parallel approach. Computers and geosciences, 2002, 28, 9:1069-1077.
2. Guo X, Lange M, Gorman G, Mitchell L, Weiland M. Developing a scalable hybrid MPI/OpenMP unstructured finite element model. Computers and Fluids, 2015, 110:227-234.
3. Micikevicius P. 3D finite difference computation on GPUs using CUDA. Proceedings of 2nd workshop on generalpurpose processing on graphics processing units. ACM, 2009:79-84.
4. Nickolls J, Buck I, Garland M, Skadron K. Scalable parallel programming with CUDA. Queue, 2008, 6, 2:40-53.
5. Abdelkhalek R, Calendra H, Coulaud O, Latu G, Roman J. Fast seismic modeling and reverse time migration on a GPU cluster. HPCS'09 Intern. Conf. on High Performance Computing and Simulation, IEEE, 2009:36-43.
6. Foltinek D, Eaton D, Mahovsky J, Moghaddam P, McGarry R. Industrial-scale reverse time migration on GPU hardware. 2009 SEG Annual Meeting. - Society of Exploration Geophysicists. 2009.
7. Bohlen T. Parallel 3-D viscoelastic finite difference seismic modelling. Computers and Geosciences, 2002, 28, 8:887-899.
8. Gropp W, Lusk E, Doss N, Skjellum A. A highperformance, portable implementation of the MPI message passing interface standard. Parallel computing, 1996, 22, 6:789-828.
9. Martin R, Komatitsch D, Blitz C, Le Goff N. Simulation of seismic wave propagation in an asteroid based upon an unstructured MPI spectral-element method: blocking and non-blocking communication strategies. Intern. Conf. on High Performance Computing for Computational Science. Springer, Berlin Heidelberg, 2008:350-363.
10. Rostrup S, De Sterck H. Parallel hyperbolic PDE simulation on clusters: Cell versus GPU. Computer Physics Communications, 2010, 181, 12:164-179.
11. Aochi H, Dupros F. MPI-OpenMP hybrid simulations using boundary integral equation and finite difference methods for earthquake
dynamics and wave propagation: Application to the 2007 Niigata Chuetsu-Oki earthquake. Procedia Computer Science, 2011, 4:1496-1505.
12. Vanderbauwhede W Benkrid K. High-performance computing using FPGAs. New York, Springer, 2013.
13. Krueger J, Donofrio D, Shalf J, Mohiyuddin M, Williams S, Oliker L, Pfreundt FJ. Hardware/ software co-design for energy-efficient seismic modeling. Proceedings of 2011 Intern. Conf. for High Performance Computing, Networking, Storage and Analysts. ACM, 2011, 73 p.
14. Petrov IB, Favorskaya AV, Sannikov AV, Kvasov IE. Grid-characteristic method using high-order interpolation on tetrahedral hierarchical meshes with a multiple time step. Mathematical Models and Computer Simulations, 2013, 5, 5:409-415.
15. Golubev VI, Petrov IB, Khokhlov NI. Numerical simulation of seismic activity by the grid-characteristic method. Computational Mathematics and Mathematical Physics, 2013, 53, 10:1523-1533.
16. Beklemysheva KA, Petrov IB, Favorskaya AV. Numerical simulation of processes in solid deformable media in the presence of dynamic contacts using the grid-characteristic method. Mathematical Models and ComputerSimuhtions, 2014, 6, 3:294-304.
17. Butenhof DR. Programming with POSIX threads. Addison-Wesley Professional, 1997, 400 p.
18. Munshi A. The opencl specification. 2009 IEEE Hot Chips 21 Symposium (HcS). IEEE, 2009, 314 p.
19. LeVeque RJ. Finite volume methods for hyperbolic problems. Cambridge university press, 2002, 558 p.
20. Carcione JM, Herman GC, ten Kroode APE. Review Article: Seismic modeling. Geophysics, 2002, 67, 4:1304-1325.
Иванов Андрей Михайлович
студент
Московский физико-технический институт Долгопрудный, 141701 Московская обл., Россия ip-e@mail.ru
Хохлов Николай Игоревич
к.ф.-м.н, с.н.с.
Московский физико-технический институт
Долгопрудный, 141701 Московская обл., Россия
k_h@inbox.ru
Петров Игорь Борисович
д.ф.-м.н, проф., член-корр. РАН и РАЕН
Московский физико-технический институт
Долгопрудный, 141701 Московская обл., Россия
НИИ системных исследований РАН
36/1, Нахимовский просп., 117218 Москва, Россия
petrov@mipt.ru
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
THE USE OF PARALLEL PROGRAMMING TECHNOLOGIES FOR MODELING SEISMIC WAVES USING GRID-CHARACTERISTIC METHOD
1Andrey M. Ivanov, ^ikolay I. Khokhlov, 12Igor B. Petrov
Moscow Institute of Physics and Technology, http://mipt.ru 141701 Dolgoprudny, Moscow region, Russian Federation
^Scientific Research Institute of System Analysis, Russian Academy of Sciences, https://www.niisi.ru 117218 Moscow, Russian Federation ip-e@mail.ru, k_h@inbox.ru, petrov@mipt.ru
Abstract. The paper deals with the problem of modeling seismic wave propagation in solid deformable bodies. The state of these bodies can be represented as a solution of hyperbolic system of equations. Grid-characteristic method is used to find the solution. This formulation of physical problem is used to demonstrate the use of technologies of parallel programming for shared memory systems using OpenMP and POSIX Threads technology. To achieve high performance, vector SSE and AVX instructions of CPU is used. In addition, GPU implementation is written using CUDA and OpenCL technologies. Results of multi-GPU parallelisation is also obtained. Use of NVIDIA GPUDirect technology that allows GPUs to make data exchanges directly through the PCI Express bus is considered. The paper demonstrates achieved acceleration and percentage of the theoretically possible performance. The paper demonstrates achieved acceleration and percentage of the theoretically possible performance. The paper describes the optimizations, which allowed us to achieve the obtained acceleration results.
The study was financially supported by RFBR as part of a research project number 15-07-01931 A.
Keywords: mathematical modeling, parallel programming, grid-characteristic method, shared memory, seismic, vectorization, GPU, CUDA, OpenCL
UDC 519.633
Bibliography — 20 references RENSIT, 2016, 8(2)185-195
Received 25.11.2016 DOI: 10.17725/rensit.2016.08.185
