Разработка алгоритмов решения систем гиперболических уравнений на графических процессорах (1187418), страница 3
Текст из файла (страница 3)
Рис. 4.2
Как видно из Рис. 2.1 результат получился таким же, как и в случае последовательного алгоритма. Тут размер сетки взят чуть больше, а именно 1000x1000. Здесь показана картина через 20 шагов.
Попробуем выяснить, насколько оптимален алгоритм, для этого, как и в случае последовательного алгоритма построим зависимость времени, приходящегося на расчет одного узла сетки, от размера сетки.
График 2.1
Как видно из График 2.1 зависимость можно считать линейной. Это придает уверенности в том, что алгоритм написан адекватно с точки зрения скорости работы и масштабируемости.
Теперь рассмотрим, какого ускорения удалось добиться.
График 2.2
На График 2.2 можно увидеть как зависит ускорение в зависимости от характерного размера сетки, в данном случае это просто квадратная сетка, со стороной, равной характерному размеру. Здесь измерялось время работы программы без инициализации. Под инициализацией понимается процесс подготовки графического устройства к вычислениям (например: копирование данных с host на deice). Алгоритм тестировался на различных CPU, которые также приведены на графике.
Из График 2.2 можно сделать вывод, что ускорение увеличивается, тем больше, чем больше размер сетки. Это можно объяснить тем, что все же некоторое время уходит на запуск вычислений на графическом устройстве, которые нельзя предотвратить.
Важно отметить тот факт, как измерялось ускорение, а точнее время работы на GPU. Если с CPU есть много стандартных, точных методов измерения времени, то для GPU важно не пропустить тот факт, что вызовы функций, которые должны выполняться на GPU, происходят асинхронно с работой CPU. Т.е. при вызове функции управление сразу возвращается к CPU. Соответственно, время работы, если не синхронизировать CPU и GPU будет считаться не верно. Поэтому во всех временных тестах перед остановкой таймера вызывается функция синхронизации.
Сразу коснемся еще одной проблемы, возникающей как раз из-за того, что крайней невыгодно часто копировать информация с/на GPU. Поэтому получается, что все вычисления надо проводить с достаточно небольшим объемом данных. По крайней мере, для практических задач нужны сетки, в которых содержится до узлов. В нашем распоряжении имеется только 2GB памяти на GPU. Задача переноса достаточно легкая, в том смысле, что каждый узел содержит мало информации. Максимальный размер сетки при решении уравнения переноса получился 10000x10000 узлов. Что, конечно нельзя сказать о решении системы гиперболических уравнений, которую мы рассмотрим далее.
-
Уравнение упругости
-
Постановка задачи
-
Уравнение упругости играет большую роль в моделировании многих физических процессов. Например, при разведке углеводородных месторождений. Т.к. на данный момент разведка, в основном, осуществляется с помощью взрывов на поверхности, которые в свою очередь создают акустические волны, которые мы и моделируем с помощью уравнения упругости.
В зависимости от условий динамического деформирования и от реологии материала в механике деформируемого твердого тела существуют различные замкнутые математические модели.
Используемая система двумерных динамических уравнений в тензорной форме имеет вид
- уравнения движения, (28)
- реологические соотношения. (29)
Здесь ρ – плотность среды, vi – компоненты вектора скорости, σij и εij – компоненты тензоров напряжения и деформаций, - ковариантная производная по k-ой координате,
– добавочная правая часть. Тензор четвертого ранга
задается в соответствии с реологией среды. В приближении малых деформаций (линейно-упругое тело) его можно записать в виде
(30)
где – символ Кронекера,
и
– постоянные Ламе.
Использовалось матричное представление данных уравнений [6]:
(31)
В этом представлении u = (v1, v2, σ11, σ12, σ22)T – вектор искомых функций, f – вектор правых частей. Явный вид матриц , а также аналитическое выражение их собственных значений и собственных векторов приведены в [3].
-
Описание схемы и ограничителя, который используется в расчетах.
В данной части работы использовалась TVD-разностная схема [3] 2-го порядка точности, описанная в параграфе 1.4. Но в данной схеме, мы будем использовать другой ограничитель.
Основной ограничитель, использованный во всех расчетах, – superbee [3], он имеет следующий вид:
(32)
-
Распараллеливание решения 2-х мерной задачи упругости с помощью технологии CUDA.
Рассмотрим систему гиперболических уравнений:
(33)
Т.к. система уравнений гиперболическая, это означает, что матрица A – невырожденная. Тогда ее можно представить в виде:
(34)
Где Ω – матрица, составленная из собственных векторов матрицы , а Λ – диагональная матрица. Тогда, домножив обе части равенства на матрицу Ω слева, получим:
(35)
Теперь обозначим . Тогда получим систему:
(36)
В нашем, двумерном случае, мы получим:
(37)
(38)
(39)
(40)
(41)
(42)
(43)
(44)
Тогда можно получить систему:
(45)
В результате можно получить явный вид матриц и
:
(46)
(47)
Приведем также явный вид матриц и
:
(48)
Где – поперечная скорость звука в среде,
– продольная скорость звука в среде.
(49)
Аналогично для матриц и
.
Как видно, данная система приведена к диагональному виду, а значит, можно решать данную систему построчно на каждом шаге, т.к. на данном шаге уравнения независимы. И решая последовательно каждый такой шаг можно сильно распараллелить. В 2-х мерном случае таких шага 2, сначала считаем, как изменилась система по оси x, потом, как изменилась система по оси y.
Здесь так же хочется заметить, что алгоритм построен таким образом, что нет необходимости на каждом шаге обмениваться данными с хостом.
Для того чтобы оттестировать программу, необходимо взять тестовый расчет. В качестве такового было взято моделирование точечного взрыва в центре квадратной сетки. Результаты распараллеленной программы совпадают с результатами последовательной.
Программа работает 1000 шагов, причем результат сохраняется каждые 10 шагов.
Сетка сделана таким образом, что граница сетки полностью поглощает возмущение. Поэтому никаких отражений не наблюдается. Идейно это реализовано, как копирование последних 2-х вычисляемых значений сетки на предыдущем шаге, для подсчета крайней точки (потому что схема основывается на 5 точках), т.е. происходит продолжение, а значит, отраженной волны не возникает.
Рис. 5.1
На Рис. 3.1 изображен результат работы программы через 200 шагов работы программы.
Рис. 5.2
На Рис. 3.2 изображен результат работы программы через 800 шагов.
На рисунках Рис. 3.2 Рис. 3.1 явно видно распространение фронта точечного взрыва.
Идея распараллеливания 2-х мерной задачи состоит в следующем. Каждый поток обрабатывает свой узел. В начале, рассмотрим алгоритм, который конечно ускоряет работу программы, но еще не использует большую часть преимуществ технологии CUDA.
В новой задаче каждый узел становится “тяжелее”, т.е. в нем появляется новая информация, новые данные, которые надо хранить. Вследствие чего, и изложенных в главе 2.2 обстоятельств, максимально возможный размер сетки уменьшается.
Для сравнения с CPU использовалась готовая реализация алгоритма результаты работы параллельной программы мы будем сравнивать именно с ним. Потому что он был оттестирован и максимально ускорен.
Как уже говорилось ранее, в начале идея работы заключалась в том, что каждый поток обрабатывает 1 узел, вот какие результаты по ускорению были получены.
График 3.1
На График 3.1 изображена зависимость ускорения, в зависимости от размера сетки. Можно сделать вывод, что ускорение либо увеличивается, при увеличении размера сетки, либо не меняется, безусловно, есть смысл сравнивать с наименьшим временем работы.
Как результат, зафиксируем, что ускорение более 10 раз, в зависимости от размера сетки. В [7] написано, что это стандартный результат для базового распараллеливания любого последовательного алгоритма, если мы не используем преимуществ GPU, а точнее, не опираемся на внутренне устройство GPU.
Но, сейчас мы считаем так, что каждый поток, обсчитывает 1 узел и не использует разделяемую память. Поэтому из логических соображений, довольное большое время уходит на получение данных из общей памяти. Т.к. количество узлов, для подсчета нового значение при схеме TVD составляет 5, т.е. надо знать 2 узла справа и 2 слева (если расчет ведется по оси x). Поэтому каждый процесс независимо будет доставать из общей памяти 5 значений.
Здесь хочется сделать ремарку, считывание из памяти устроено таким образом, что 1 значение никогда не считывается, а считывается сразу несколько подряд лежащих значений, это называется bulk. Поэтому всегда быстрее считывать данные именно по порядку, а не “против шерсти” или беспорядочно. Хотя для расчета вдоль оси X, это не играет большой роли, то при подсчете по оси Y получается такая ситуация, что каждое такое значение лежит на довольно большом расстоянии друг от друга, а значит, нет шанса считать эти значения подряд за один раз, а приходится 5 раз лазить за данными.
Учитывая этот факт, становится понятно, что количество узлов, обрабатываемых процессором, не может быть равно одному. Поэтому было проведение исследование, для выяснения того, какое лучшее число узлов.
График 3.2 Сетка 128x128, 1000 шагов