Software & Systems

В процессе решения задач механики деформируемого твердого тела при конечных деформациях для сжимаемых и несжимаемых материалов методом конечных элементов получаются системы нелинейных уравнений. Их решение методом Ньютона–Рафсона сводится к решению на каждой итерации системы линейных алгебраических уравнений (СЛАУ) с разреженной матрицей (размерность матрицы зависит от размерности задачи и того, насколько мелкая сетка используется в методе конечных элементов) с нерегулярной структурой. Решение задачи тем точнее, чем мельче сетка, поэтому актуальна возможность решения системы уравнений максимально возможной размерности при заданных ресурсах ЭВМ [1–3].

В работе реализован программный модуль, объединивший различные алгоритмы и методы: прямые и итерационные решатели для симметричных и несимметричных матриц системы, различные предобуславливатели в итерационных методах и способы хранения матрицы в памяти, параллельные вычисления с использованием технологий OpenMP и CUDA [4–7]. В программном модуле реализован метод решения задач для несжимаемых материалов на основе алгоритма Узавы [8–12].

Для программного модуля разработан и реализован алгоритм оптимального выбора решателя в зависимости от механической постановки задачи, ее размерности и возможностей компьютера [13]. При желании пользователь может сам ограничивать некоторые возможности выбора и задавать параметры, влияющие на выбор решателя, или даже указать явно, какой решатель он хочет использовать. По сути программный комплекс является некоторой оболочкой над отдельными решателями, которая принимает от программы, выполняющей механическую задачу (будем называть ее ядром), матрицу системы, правую часть и при желании некоторые параметры настройки, а затем в рамках содержащегося в ней алгоритма определяет, какой именно решатель необходимо запускать, настраивает его и приводит матрицу к соответствующему виду (разные решатели могут иметь разные оптимальные форматы хранения для матриц). После окончания работы решателя результат решения возвращается в ядро для дальнейшей обработки. На рисунке 1 схематично показан принцип работы модуля и его взаимодействия с ядром.

Основные учитываемые критерии и их влияние на выбор решателя:

1)    если матрица малой размерности (критерий основывается на доступных ресурсах используемого компьютера), выбирается решатель, использующий прямой метод;

2)    если задача имеет симметричную матрицу, в качестве итерационного метода целесообразно использовать метод сопряженных градиентов (CG), в противном случае – бисопряженных градиентов (BiCG) или обобщенных минимальных невязок (GMRes);

3)    используется наиболее эффективный (быстрый, но необязательно самый «сильный») предобуславливатель, имеющийся для выбранного решателя;

4)    если компьютер позволяет использовать технологию CUDA, выбирается решатель, использующий эту технологию;

5)    если в задаче используется несжимаемый материал, выбирается метод на основе алгоритма Узавы.

Проведен ряд численных экспериментов, подтверждающих обоснованность приведенных выше критериев.

Если по какой-то причине после выбора решателя задача не была решена (например, если выбранным методом и с выбранным предобуславливателем итерационный процесс не сошелся), может быть запущен другой вариант (например, с более «сильным» предобуславливателем).

Программный модуль уже был использован в программном комплексе прочностного анализа CAE FIDESYS [14]. Пользовательский интерфейс для задания настроек программного модуля создан так, что некоторые параметры (например, указывающий на решение задачи с несжимаемым материалом или на симметричность задачи) скрыты от пользователя и определяются программно в зависимости от постановки задачи или параметров, заданных пользователем для материала.

Выбор решателя 1. Всегда, если есть возможность, стоит использовать прямые методы. Их преимущество в том, что они дают точное решение (итерационные методы будут давать приближенное), конечно, с учетом машинной точности. Другим преимуществом прямых методов является то, что они находят решение за один проход, поэтому очень быстры. Но у них есть и серьезный недостаток: требуется большой объем памяти для проведения численной факторизации, поэтому при больших размерах матрицы их использование становится невозможным. Если, помимо прочего, известно, что матрица системы симметрична, то стоит использовать это для оптимизации (например, можно экономить память, храня только верхнюю часть матрицы) [4, 5].

Выбор решателя 2. В большинстве реальных 3D-задач матрица имеет слишком большую размерность и оперативной памяти не хватает для того, чтобы можно было пользоваться прямыми методами. В таких случаях используются итерационные методы (проекционные методы на сопряженных подпространствах пространства Крылова, основанные на минимизации невязки).

Здесь симметричность может служить уже не просто для оптимизации. Существуют специальные методы только для симметричных СЛАУ, которые являются более стабильными и эффективными, нежели методы для несимметричных систем. В симметричном случае хорошо зарекомендовал себя метод сопряженных градиентов (Conjugate Gradient). Для несимметричных задач можно использовать его модификацию (Biconju­gate Gradient Stabilized) или модификацию метода обобщенных минимальных невязок (Flexible Gene­ralized Minimal Residual) [4–7].

Выбор решателя 3. При сложной геометрии задачи и использовании неструктурированных расчетных сеток портрет матрицы может быть нерегулярным (портрет – множество пар индексов, соответствующих ненулевым элементам), а число ее обусловленности – очень большим (от этой величины зависит, насколько хорошо сходится решение задачи итерационными методами; чем больше эта величина, тем решение сходится хуже).

Для уменьшения числа обусловленности в программном модуле реализованы различные возможности использования предобуславливателей [4–7]:

–      без предобуславливателя;

–      предобуславливатель Якоби (диагональный);

–      SSOR (symmetric successive overrelaxation – симметричный метод последовательной верхней релаксации);

–      AINV – параметризуемый, поэтому реализовано несколько видов под разные наборы параметров;

–      неполное LU-разложение и его разновидности: ILU(0), ILU(2), ILUT;

–      AMG (алгебраический многосеточный).

Стоит отметить, что некоторые из них подходят только для симметричных случаев, так как использование такого предобуславливателя может нарушить симметрию СЛАУ. Некоторые предобуславливатели реализованы только для CPU, а другие эффективны только для GPU.

Обратим также внимание на то, что использование предобуславливателя вынуждает нас выполнять некоторые дополнительные действия, связанные с созданием предобуславливателя и домножением на него. Чтобы убедиться в оправданности этих действий и в том, что уменьшение количества итераций перекрывает все дополнительные затраты на использование предобуславливателя, были проведены некоторые численные эксперименты. Как правило, чем сложнее в построении предобуславливатель, тем больше времени требуют эти дополнительные действия. Но есть и обратная сторона: чем сложнее предобуславливатель, тем сильнее он уменьшает число обусловленности матрицы, а значит, дает лучшую сходимость. Необходимо искать компромисc.

Выбор решателя 4. Итерационные методы хорошо поддаются распараллеливанию. Поэтому, если позволяет архитектура компьютера, очень важно использовать такую возможность. Многопоточность на центральном процессоре (CPU) есть уже практически на любом современном компьютере. Но если компьютер позволяет осуществлять распараллеливание счета с использованием видеокарты (GPU), это может давать гораздо более заметный прирост производительности и не стоит упускать такую возможность. Однако это касается не всех видеокарт. Одними из наиболее популярных являются видеокарты NVIDIA, в современных моделях которых имеется поддержка технологии CUDA. Именно эта технология и позволяет перенести решение СЛАУ с CPU на GPU [13].

На рисунке 2 сравниваются результаты расчета на центральном процессоре (Host) и на графическом устройстве (Device) на двух различных персональных компьютерах (AMD Athlon 64 X2 Dual Core с графическим устройством NVIDIA GF GTX 260 и Intel Pentium Dual Core E6500 с графическим устройством NVIDIA GF GTX 295). Расчет делался для задачи теории упругости с симметричной матрицей системы линейных алгебраических уравнений, содержащей 300 тыс. строк; расчеты производились методом сопряженных градиентов без предобуславливателя и с предобуславливателем Якоби (на рисунке обозначено как Jacobi). Как видно из рисунка, использование графического устройства дает заметно большее преимущество, чем увеличение количества процессоров. Предобуславливание тоже дает очень заметный прирост. А если их совмещать, время расчета значительно уменьшается (в 40–75 раз).

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

Известно несколько вариантов алгоритма Узавы, основанных на различных релаксационных методах решения СЛАУ. В данном программном модуле реализованы варианты, основанные на методе наискорейшего градиентного спуска, методе сопряженных градиентов (двух- и трехслойная схемы) и трехслойном методе сопряженных невязок. Расчетные формулы для этих вариантов метода Узавы записываются по аналогии с расчетными формулами соответствующих итерационных методов [8–12].

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

Сравним полученные результаты для четырех матриц различной размерности:

1)    30 402 строки, из них 20 402 приходятся на главный блок (матрица А);

2)    120 802 строки, из них 80 802 приходятся на главный блок;

3)    246 534 строки, из них 164 738 приходятся на главный блок;

4)    481 602 строки, из них 321 602 приходятся на главный блок.

Критерием окончания расчета было уменьшение нормы невязки в 100 раз по сравнению с нормой начальной невязки. Для решения СЛАУ с матрицей А на каждой итерации метода Узавы использовались прямые методы.

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

Отметим, что количество итераций алгоритма Узавы почти не зависит от размерности матрицы, но, тем не менее, как можно видеть из графика, время расчета увеличивается с размером матрицы (матрицы на графике упорядочены по возрастанию их размерности). Это связано с ростом времени, потраченного на одну итерацию алгоритма и решение СЛАУ с матрицей A, размерность которой тоже растет. Также можно видеть, что наиболее эффективным и стабильным является метод Узавы, основанный на методе сопряженных невязок (CRes).

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

Литература

1.     Левин В.А., Калинин В.В., Зингерман К.М., Верши- нин А.В. Развитие дефектов при конечных деформациях. Компьютерное и физическое моделирование. М.: Физматлит, 2007. 392 с.

2.     Zienkiewicz O.C., Taylor R.L. The finite element method. The basis. Oxford: Butterworth-Heineman, 2000, vol. 1, 691 p.

3.     Levin V.A., Vershinin A.V. Non-stationary plane problem of thesuccessive origination of stress concentrators in a loaded bo­dy. Finite deformations and their superposition. Communications in Numerical Methods in Engineering, 2008, vol. 24, pp. 2229–2239.

4.     Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 2003. 632 с.

5.     Богачев К.Ю. Методы решения линейных систем и нахождения собственных значений. М.: Изд-во мех.-мат. ф-та МГУ, 1998.

6.     Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. Н.: Изд-во НГТУ, 2000.

7.     Saad Y. Iterative Methods for Sparse Linear Systems, Second edition, SIAM, Philadelphia, 2003.

8.     Быченков Ю.В., Чижонков Е.В. Итерационные методы решения седловых задач. М.: БИНОМ. Лаборатория знаний, 2010. 349 с.

9.     Чижонков Е.В. Релаксационные методы решения седловых задач. М.: Изд-во ИВМ РАН, 2002. 238 с.

10.  Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970.

11.  Быченков Ю.В. Оптимизация предобусловленных методов для седловых задач // Докл. Акад. наук. 2002. Т. 384. № 4. С. 439–441.

12.  Степин Н.Е., Левин В.А., Зингерман К.М., Верши- нин А.В. Сравнительный анализ различных вариантов алгоритма Узавы в задачах упругости для несжимаемых материалов // Вестн. ТвГУ: Сер. Прикладная математика. № 3 (26). С. 29–34.

13.  Левин В.А., Вершинин А.В., Траченко А.В., Прокопенко А.С., Степин Н.Е. Некоторые результаты использования технологии CUDA при решении СЛАУ для задач прочности при перераспределении конечных деформаций / В кн.: Современные проблемы математики, механики, информатики: матер. 10-й Междунар. конф. Тула, 2009.

14.  Levin V.A., Zingerman K.M., Vershinin A.V., Nikifo- rov I.V. CAE FIDESYS for strength analysis at large strains and their redistribution. 10th Word Congress on Computational Mechanics, 8–13 July 2012, Sao Paulo, Brazil, 2012, 19579, p. 323.

References

1.  Levin V.A., Kalinin V.V., Zingerman K.M., Vershinin A.V.  Razvitie defektov pri konechnykh deformatsiyakh.
Kompyuternoe i fizicheskoe modelirovanie [Defects  Development at Finite Deformations. Computational and  Physical Mod-eling]. Moscow, Fyzmatlit Publ., 2007, 392 p.

2.  Zienkiewicz O.C., Taylor R.L.  The finite element method. Vol.  1.  The basis. Oxford,  Butterworth-Heineman  Publ.,
2000, 691 p.

3.  Levin V.A., Vershinin A.V. Non-stationary plane problem of the successive origination of stress concentrators in a
loaded body. Finite deformations and their superposition.  Communications in Numerical Methods in Engineering.  2008,
no. 24, pp. 2229–2239.

4.  Bakhvalov N.S.,  Zhidkov N.P., Kobelkov G.M.  Chislennye metody  [Numerical  Methods].  Moscow,  Science  Publ.,
2003, 632 p.

5.  Bogachov K.U. Methody resheniya lineynykh system i nakhozhdeniya sobstvennykh znacheniy [Linear System Solv-ing Methods and Eigenvalues Finding]. Moscow, MSU Mechanic and Mathematic Department Publ., 1998.

6.  Balandin M.Yu., Shurina E.P.  Metody resheniya SLAU bolshoy razmernosti  [Solving  Methods of  Large Dimension
Linear Algebraic Equation Systems]. Novosibirsk, NNSTU Publ., 2000.

7.  Saad Y. Iterative methods for sparse linear systems. 2nd edition, SIAM, Philadelphia, 2003.

8.  Bychenkov  Yu.V., Chizhonkov E.V.  Iteratsionnye metody resheniya sedlovykh zadach  [Iterative  Methods for  Solv-ing Saddle Point Problems]. Moscow, BINOM Knowledge Laboratory Publ., 2010, 349 p.

9.  Chizhonkov E.V. Relaksatsionnye metody resheniya sedlovykh zadach [Relaxation Methods for Solving Saddle Point
Problems]. Moscow, 2002.

10.  Ladyzhenskaja O.A. Matematicheskie voprosy dinamiki vyazkoy neszhimaemoy zhidkosti  [The  Mathematical Theory
of Viscous Incompressible Fluid]. Moscow, Science Publ., 1970.

11.  Bychenkov Yu.V. Preconditioner methods optimization for saddle point problems.  Dokl. Akad. Nauk  [Reports of the
Academy of Sciences]. 2002, vol. 384, no. 4. pp. 439–441 (in Russ.).

12.  Stepin N.E., Levin V.A., Zingerman K.M., Vershinin A.V. Comparative analysis of Uzava algorithm various options
for elasticity problems of incompressible materials.  Vestn. TvGU: Ser. Prikladnaya matematika  [Bulletin of the Tver State
Univ. Applied Mathematics Series], no. 3 (26), pp. 29–34 (in Russ.).

13.  Levin V.A., Vershinin A.V., Trachenko A.V., Prokopenko A.S., Stepin N.E. Some results of using CUDA technolo-gy for linear algebraic equation systems solving for problems of strength in the redistribution of finite deformations.  Mater.

10 Mezhdunar. konf. “Sovremennye problemy matematiki, mekhaniki, informatiki”  [Proc.  10th  Int.  Conf.  Modern Problems
of Mathematics, Mechanics, Informatics]. Tula, 2009.

14.  Levin V.A., Zingerman K.M.,  Vershinin A.V., Nikiforov I.V.  CAE FIDESYS for strength analysis at large strains
and their redistribution. Proc. 10th Word Congress on Computational Mechanics. 8–13 July 2012, Sao Paulo, Brazil, 19579,
p. 323.