Применение параллельного алгоритма BiCGStab для решения двухмерной упругопластической задачи

Применение параллельного алгоритма BiCGStab для решения двухмерной упругопластической задачи[1]

, ,

Россия, Екатеринбург

Введение

Упругопластическая задача с большими пластическими деформациями физически и геометрически существенно нелинейная и требует большого количества времени для её решения методом конечных элементов (МКЭ) на персональном компьютере. Для сокращения времени расчётов необходимо использовать параллельные вычисления, в частности, на кластерных системах.

Реализация МКЭ в упругопластических задачах осуществляется в условиях пошагового нагружения и на каждом таком шаге состоит из трёх основных этапов [1]:

1) расчёт локальных матриц жёсткости для конечных элементов и формирование глобальной матрицы жёсткости и вектора правой части системы линейных алгебраических уравнений (СЛАУ)

(1)

относительно искомого вектора обобщённой скорости в узлах конечно-элементной сетки;

2) решение СЛАУ (1);

3) вычисление напряжённо-деформированного состояния конечных элементов в конце

шага нагружения.

На каждом шаге нагружения этап 1 выполняется один раз, а этапы 2 и 3 — от десяти до пятнадцати раз для удовлетворения итерационно с приемлемой точностью условию пластичности в конце шага нагружения. При этом матрица жёсткости остаётся постоянной, а изменяется только правая часть системы уравнений. Если этапы 1 и 3 легко распараллеливаются, то алгоритмизация распределённого параллельного решения СЛАУ является сложной задачей.

Матрица имеет ленточный вид с небольшой относительной шириной ленты. Ширина ленты зависит как от типа аппроксимирующей расчётную область сетки, так и от способа нумерации её узлов и является единственной основной характеристикой матрицы СЛАУ. Для рассматриваемого алгоритма решения системы уравнений в данной работе, с целью упрощения анализа результатов, используются только регулярные сетки с одинаковым числом ячеек вдоль каждой из сторон. Для таких сеток отношение полуширины ленты к размерности вектора (к количеству переменных в задаче) для двумерных равно 0,01, 0,005 и 0,0033 соответственно для количества разбиений конечноэлементной сетки =100, 200 и 300

Ранее была проведена работа по исследованию применимости итерационных методов решения СЛАУ для двумерной упругопластической задачи с ленточной матрицей [2]. В отличие от работы [2], в данной работе рассматривается разряжённые схемы хранения матрицы жёсткости и используется стабилизированный метод бисопряжённых градиентов BiCGStab [3], который позволяет решать СЛАУ с несимметричной матрицей.

Решение СЛАУ

При решении задачи использовали две разряжённые схемы хранения матрицы жёсткости. Для построения глобальной матрицы жёсткости использовали координатную схему хранения матрицы [4]. В данном случае вместе с элементом матрицы хранятся его координаты в матрице. Эти тройки чисел хранятся в самобалансирующемся дереве поиска [5], с заданной сортировкой по координатам в матрице. Такое хранение даёт возможность эффективно добавлять элементы в произвольное место матрицы при её построении.

Для проведения вычислений с матрицей использовали сжатую по строкам форму хранения матрицы (CRS) [4]. Для этого формата строятся три массива:

·  - в котором хранятся значения элементов матрицы. Размер массива равен количеству ненулевых элементов матрицы жесткости.

·  - этот массив содержит номера элементов в массиве, являющимися первыми ненулевыми элементами в строках матрицы жёсткости. Последним элементом массива является количество ненулевых элементов матрицы жёсткости

·  - в котором хранятся номера столбцов ненулевых элементов матрицы жесткости. Количество элементов массива совпадает с количеством элементов массива .

Форма CRS удобна и эффективна для проведения вычислений с матрицами, особенно операций произведения матрицы на вектор, составляющих большую часть итерационных методов. С другой стороны форма CRS неэффективна при добавлении нового элемента внутрь матрицы, так как при этом требуется сдвигать все элементы в массивах и и перестраивать массив .

Предобуславливатели

Для повышения скорости сходимости решения СЛАУ (1) итерационными методами применяют предобуславливатель [4]. Он представляет собой квадратную невырожденную матрицу. Размерность предобуславливателя совпадает с размерностью матрицы системы. В работе применили диагональный предобуславливатель, состоящий из диагональных элементов матрицы . Он позволил сократить количество итераций приблизительно 2 раза по сравнению с непредобусловленным вариантом алгоритма.

Стабилизированный метод бисопряжённых градиентов

Стабилизированный метод бисопряжённых градиентов является модификацией метода сопряжённых градиентов для несимметричных матриц. При решении СЛАУ использовали модификацию метода BiCGStab для машин с распределённой памятью, описанную в работе [6]. Данная модификация позволяет проводить межпроцессорный обмен, требующийся для вычисления скалярных произведений, одновременно с применением предобуславливателей.

Вычислительный эксперимент

Вычислительные эксперименты проводились на кластере umt Института математики и механики УрО РАН. Каждый из 208 вычислительных узлов оборудован двумя 4-х ядерными процессорами Intel Xeon E5450 и 16 ГБ оперативной памяти. Узлы соединены между собой Infiniband для передачи сообщений и Gigabit Ethernet для ввода-вывода. Для компиляции использовался компилятор Intel версии 1Для передачи сообщений использовался OpenMPI.

Анализ эффективности применения параллельных вычислений был произведён на примере решения МКЭ двумерной задачи сжатия цилиндра плоскими плитами из упругопластического изотропного и изотропно-упрочняемого материала. Решение основывается на принципе виртуальной мощности в скоростной форме [1]:

(2)

со следующими определяющими соотношениями, полученными в работе [7]:

,

, , , ,

, - условие текучести Мизеса (3)

.

Здесь - тензор напряжений Коши; - плотность поверхностных сил; - промежуток времени для шага приращения нагрузки; - вариация кинематически допустимых полей скоростей; - набла-оператор; , - объем и поверхность тела соответственно; , - элементы объёма и площади поверхности цилиндра соответственно; , - коэффициенты Ламе; - единичный тензор; точкой и двумя точками обозначено соответственно скалярное и двойное скалярное произведение тензоров; точкой сверху обозначена полная производная по времени; - градиент скорости перемещений; - тензор скоростей деформаций; - девиатор тензора напряжений; - среднее нормальное напряжение, - объёмный модуль упругости; - относительное изменение индивидуального объёма бесконечно малой частицы среды; - напряжение текучести; - числовой коэффициент, принимающий значения при или при , и при , ; - параметр упрочнения, в силу малости упругих деформаций .

На контакте с плитами приняли условие прилипания металла к плитам. Нагрузку в виде перемещения плиты прикладывали малыми шагами . Шаг выбирали так, чтобы отношение ( - высота параллелепипеда или цилиндра) не превышало предела упругости по деформации, в нашем случае 0,002, что обеспечивало устойчивость вычислительной процедуры. Величину относительного сжатия параллелепипеда и цилиндра приняли равной 0,5.

На каждом шаге нагрузки равенство (2) с помощью конечно-элементной аппроксимации сводится к СЛАУ (1).

При распараллеливании формирования глобальной матрицы жёсткости конечные элементы распределялись слоями между процессорами. При этом операция объединения матрицы жёсткости осуществлялась только с соседними процессорами. При таком распараллеливании количество процессоров не может превышать количество слоёв конечных элементов.

Параметры матрицы для вычислительного эксперимента представлены в таблице.

Таблица. Параметры сеток в вычислительных экспериментах в задаче сжатия цилиндра

100	205	20402	0,0100
200	405	80802	0,0050
300	605	181202	0,0033

Зависимости ускорения от количества процессоров при решении СЛАУ показаны на рисунке. Под ускорением понимается величина , где и - время, затраченное на вычисление на и 1 процессорах соответственно.

Рисунок. Зависимость ускорения от числа процессоров при решении СЛАУ в двухмерной упругопластической задаче.

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

Выводы

Метод BiCGStab показывает хорошую производительность при решении двухмерной упругопластической задачи. На сетке при использовании 32 процессоров достигается десятикратное ускорение по сравнению с решением задачи на одном процессоре. Снижение затрат по памяти при использовании разряжённых схем хранения позволят использовать данный метод при решении задачи на графических процессорах.

Литература

1.  Поздеев A.А., , Няшин упруго-пластические деформации. – М.: Наука, 1986. – 232 с.

2.  , , Коновалов вычислительных экспериментов решения двумерной упругопластической задачи итерационными методами на кластерной системе. // Физико-химическая кинетика в газовой динамике. 2010. Том 9. http://www. chemphys. *****/pdf/023.pdf

3.  H. A. van der Vorst, BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. put.31–644. doi:10.1137/0913035.

4.  Saad Y. Iterative methods for sparse linear systems, SIAM 2003

5.  Кормен, Лейзерсон, Ривест, Клиффорд Штайн «Алгоритмы: построение и анализ», 2-е здание 2008 М. . Д. Вильямс», 1296 с.

6.  Boris Krasnopolsky, The reordered BiCGStab method for distributed memory computer systems. Procedia Computer Science, Volume 1, Issue 1, ICCS 2010, May 2010, p. 213-218

7.  Коновалов соотношения для упругопластической среды при больших пластических деформациях // Известия РАН. Механика твердого тела. 1997. № 5. С. 139-147

[1] Работа выполнена по программе Президиума РАН №15, проект 12-П-1-1025