В настоящее время активно используются конечно-элементные программные комплексы (КЭПК) для моделирования деформации твердых тел методом конечных элементов (МКЭ). Применение МКЭ порождает системы линейных алгебраических уравнений (СЛАУ) большой размерности, которые необходимо решать в процессе расчёта. Матрицы таких систем, называемые матрицами жёсткости, обладают рядом свойств, зависящих от структуры и количества узлов конечно-элементной сетки, а также характера задачи. Эти матрицы в значительной степени разрежены, обладают симметричным портретом, а также блочной структурой.
Для одновременного применения сеток с большим количеством узлов и сложных моделей материалов необходимо использовать специальные подходы, позволяющие сократить время решения задачи до приемлемого уровня. Как правило, для этого используется целый ряд приёмов, включая применение современных математических методов решения СЛАУ, предобуславливателей, написания параллельного и векторизованного кода, а также задействования вычислительных ускорителей. Наиболее затратные в вычислительном плане задачи решают с использованием специальных вычислительных устройств коллективного доступа, таких как мощные серверы и суперкомпьютеры.
Современные КЭПК оснащаются набором реализаций методов решения СЛАУ. Широкая номенклатура реализаций объясняется тем, что каждая из них подходит для использования только с определённым классом задач или для определённой вычислительной системы. Следуя англоязычной терминологии, набор программ, выполняющий наиболее требовательные к производительности операции в ходе решения задачи, в данной работе будем называть «решателем» (англ. «Solver», иногда в русскоязычной литературе встречается транслитерация «Солвер»). Разработка эффективного решателя для каждой из новых вычислительных систем требует значительных усилий по проектированию архитектуры, созданию эффективного программного кода, а также проведения вычислительных экспериментов. Производители конечно-элементного программного обеспечения порой стремятся скрыть детали реализации собственного продукта от конкурентов, тем самым делая результаты внутренних исследований недоступными для научного сообщества. Фрагментарная информация о реализации методов решения СЛАУ содержится в научных публикациях и технических докладах разработчиков, озвучивается при выступлениях на конференциях или публикуется в руководствах пользователя КЭПК.
Данная статья обобщает и систематизирует подобную информацию, позволяя читателю оценить состояние проблемы, а также выделяет наиболее часто используемые подходы к решению СЛАУ в КЭПК.
Решатели в конечно-элементных программных комплексах
Можно выделить несколько характерных особенностей, которые позволяют классифицировать решатели. К ним относится применяемый метод решения СЛАУ, вычислительная система, а также структура данных для представления матрицы жёсткости.
Выделяют прямые и итерационные методы решения СЛАУ [2]. Прямые методы используют метод Гаусса, выраженный через факторизацию матрицы коэффициентов СЛАУ, а итерационные методы строят ряд приближений к решению системы.
При использовании точной арифметики прямые методы дают точное решение СЛАУ. Компьютерное представление чисел с плавающей точкой порождает зависящие от значений элементов матрицы ошибки, для уменьшения которых применяются специальные арифметические методики. Итерационные методы дают приближение к решению даже в точной арифметике, при этом скорость сходимости ряда приближений сильно зависит от числа обусловленности матрицы СЛАУ. Для ускорения сходимости методов часто используют предобуславливание [32], в ходе которого явно или неявно строится приближение матрицы, обратной к матрице коэффициентов СЛАУ.
Методы решения СЛАУ и предобуславливания
В данном разделе приводятся встречающиеся в тексте статьи методы решения СЛАУ и предобуславливания, а также даны ссылки на источники, описывающие соответствующий метод. Для некоторых методов приводятся принятые сокращения. Для того чтобы предотвратить чрезмерное разрастание списка литературы, ссылки на оригинальные статьи, описывающие методы, по возможности не приводятся. Получить более подробную информацию о методах можно с помощью приводных библиографических списков.
Прямые методы решения СЛАУ [14]: Разложение Холецкого, LDLT-разложение и LU-разложение.
Итерационные методы решения СЛАУ [32, 11]:
● SSOR (англ. «Symmetric successive over-relaxation») – метод симметричной последовательной верхней релаксации.
● CR (англ. «Conjugate Resudual») – метод сопряжённых невязок.
● CG (англ. «Conjugate Gradient») – метод сопряжённых градиентов.
● CGS (англ. «Conjugate Gradient Squared») – квадратичный метод сопряжённых градиентов.
● BiCG (англ. BiConjugate Gradient) – Метод бисопряжённых градиентов.
● BiCGStab (англ. «BiConjugate Gradient Stabilized») – стабилизированный метод бисопряжённых градиентов.
● GMRES (англ. «Generalized Minimal Residual») – обобщённый метод минимальных невязок.
● FGMRES (англ. «Flexible GMRES») – разновидность обобщённого метода минимальных невязок.
● TFQMR (англ. «Transpose-Free Quasi-Minimal Residual») – метод квазиминимальных невязок.
Методы предобуславливания, основанные на неполных разложениях [32, глава 10]:
● ILU (англ. «Incomplete LU») – неполное LU-разложение по общей схеме.
● ILU(k) – неполное LU-разложение с контролем заполнения.
● ILU(0) (англ. «Zero Fill-in ILU») – неполное LU-разложение без заполнения.
● ILUT (англ. «Incomplete LU with treshold») – неполное LU-разложение с контролем порогового значения.
● IC – неполное разложение Холецкого.
Библиотеки для решения СЛАУ
1. PARDISO
Используемые алгоритмы, а также многие детали реализации данной библиотеки описаны в научной литературе и руководстве пользователя [33]. Библиотека использует различные методы в зависимости от характеристики системы. В частности, для симметричных положительно определённых матриц используется разложение Холецкого, для симметричных матриц общего вида – LDLT-разложение, для симметричных матриц – LU-разложение. Для всех матриц выполняется сложная предобработка, включающая перестановки и масштабирование элементов, при этом возможно использование структурной симметрии несимметричной матрицы.
Помимо версии библиотеки для систем с общей памятью существует реализация для систем с распределённой памятью, основанная на MPI.
В 2006 году компания Intel начала разрабатывать собственный вариант данной библиотеки, называемой Intel PARDISO или Intel MKL PARDISO [19]. Варианты библиотеки несколько отличаются функциональными возможностями.
2. BCSLIB, BCSLIB-EXT и MUMPS
BCSLIB является коммерческой разработкой корпорации Boeing [1]. Не следует путать данную библиотеку с другой библиотекой со сходным названием, BCSLIB-EXT. Библиотека BCSLIB-EXT разрабатывается компанией Access Analytics International и имеет реализацию для графических ускорителей со сходным интерфейсом, называемую BCSLIB-GPU [27].
BCSLIB-EXT включает в себя реализации прямых методов, основанных на LU-разложении, разложении Холецкого, LDLT-разложении, алгоритмы предобработки матрицы [15].
Библиотека MUMPS предназначена для решения СЛАУ на параллельных системах с общей и распределённой памятью [22]. Используются LU- и LDLT-разложения вместе со сложной фазой предобработки, описанной в документации. Библиотека имеет возможность использовать внешнюю память.
Конечно-элементные программные комплексы
1. ANSYS
Конечно-элементный программный комплекс ANSYS является одним из старейших существующих конечно-элементных программных комплексов общего назначения. В поставку системы ANSYS входят несколько прямых и итерационных решателей. Используемые алгоритмы частично раскрываются в поставляемом вместе с документацией к комплексу справочнике по теоретическим основам [38] (англ. «Theory Reference for the Mechanical APDL and Mechanical Applications»).
Функциональность прямых решателей типична для библиотечных реализаций, основанных на факторизации матрицы: используется LU-разложение или разложение Холецкого с предобработкой. Разработчики указывают в качестве источника использованных алгоритмов монографию [28] и раскрывают конкретные алгоритмы. Первые реализации прямого решателя были разработаны внутри компании, однако с начала 2000-х годов используется библиотека BCSLIB. Решатели работают на системах с общей и распределённой памятью, поддерживают ускорители.
Первым итерационным решателем, поставляемым конечным пользователям, был решатель, обозначаемый авторами ANSYS, как JCG (a Jacobi preconditioned conjugate gradient solver) [29]. Руководство пользователя рекомендует использовать этот решатель для задач с маленькими числами обусловленности (англ. «well-conditioned»). По всей видимости, этот решатель сочетает в себе итерационный метод крыловского типа и предобуславливатель Якоби. Для симметричных матриц существует реализация решателя для нескольких графических ускорителей [12], причём ускоритель используется только для выполнения умножения матрицы на вектор.
Кроме решателя JCG вместе с ANSYS поставляется решатель PCG, также ранее называемый ANSYS PowerSolver. Авторы программы расшифровывают эту аббревиатуру как «Preconditioned Conjugate Gradient». Данный решатель лицензируется ANSYS другой корпорацией: Computational Applications and System Integration, Inc. Решатель PCG работает только с симметричными матрицами, однако позволяет решать задачи с большими числами обусловленности. Для решения СЛАУ с матрицами наиболее общего вида вместе с ANSYS поставляется решатель ICCG (Incomplete Cholesky Conjugate Gradient). Он позволяет решать СЛАУ с симметричными и несимметричными матрицами. Детали реализации данного решателя принципиально не публикуются [20]. Авторы работы [29] полагают, что в данном решателе задействуется несколько различных итерационных методов с общим обозначением. В частности, кроме метода CG с IC-предобуславливателем они сообщают об использовании методов GMRES, BiCG (нестабилизированная версия) и ILU-предобуславливатель.
2. NASTRAN
NASTRAN является одним из старейших программных комплексов конечно-элементного анализа. Существует множество версий данного конечно-элементного комлплекса, при этом вместе с различными версиями поставляются различные решатели.
Помимо оригинальной версии MSC NASTRAN, известны Autodesk NASTRAN (бывший NEi NASTRAN) и NX NASTRAN (являющийся частью Siemens PLM Software) .
Используемые в комплексе MSC NASTRAN методы описаны в руководстве пользователя [21]. Прямые решатели основаны на LU- и LDLT-разложениях и поддерживают графические ускорители, а также распределённые системы. Кроме этого, комплекс включает два итерационных решателя, основанных на методе CG, при этом один из них использует матрицы в явном виде, а второй реализует поэлементный подход.
Используемые в NX NASTRAN методы решения СЛАУ описаны в руководстве пользователя [25]. Набор прямых методов включает в себя LU-разложение, разложение Холецкого, LDLT-разложение с предобработкой для уменьшения заполненности факторизации. Используются методы CG, CR и BiCG. Для предобуславливания используется метод Якоби, IC. Методы сгруппированы в три группы реализаций, называемых в руководстве [24] решателями. Первый из них объединяет прямые методы, второй использует многосеточный метод, третий является коммерческой реализацией итерационного метода, называемый авторами «3D-iterative solver».
Вместе с Autodesk NASTRAN поставляется несколько решателей [23]: PCGLSS Solver, VIS Solver, PSS Solver и VSS Solver. PCGLSS Solver использует прямой или итерационный метод решения СЛАУ в зависимости от задачи. При решении итерационным методом используются специально разработанный предобуславливатель, использующий информацию о типе конечного элемента. В качестве итерационного метода выступает метод сопряжённых градиентов.
VIS Solver может работать с матрицами самого общего вида, в том числе и не являющимися положительно определёнными, а также имеющими высокое число обусловленности («ill-conditioned»). Решатели PSS Solver и VSS Solver используют прямые методы решения СЛАУ.
3. Autodesk simulation
Помимо встроенного в САПР компании Autodesk Nastran существует ещё один комплекс конечно-элементного анализа, принадлежащий Autodesk. Информация об используемых методах решения СЛАУ содержится в руководстве пользователя [10]. Комплекс включает четыре решателя, включая три прямых и один итерационный. К прямым решателям относятся «Banded Solver», «Skyline Solver» и «Sparse Solver». Первые два решателя используют внутренний алгоритм для оптимизации ширины ленты матрицы. Третий решатель представляет собой интерфейс для внешних библиотек решения СЛАУ, в частности BCSLIB-EXT и PVSS [26], разрабатываемый SolverSoft. Последний использует оптимизацию ширины ленты, работает с положительно определёнными матрицами и использует факторизацию. Итерационный решатель может использовать стабилизированный и обычный метод бисопряжённых градиентов, а также обобщённый метод минимальных невязок. В качестве предобуславливателей можно использовать метод релаксаций (SSOR), неполное LU-разложение с нулевым заполнением (ILU(0)) и диагональное масштабирование.
4. Abaqus FEA
Abaqus имеет два решателя, не относящихся к частотному анализу [7]. Один из решателей использует прямой метод решения СЛАУ, а другой – итерационный.
Некоторые детали реализации можно извлечь из руководства пользователя [5]. Прямой решатель («direct linear equation solver») использует прямой, разреженный вариант метода Гаусса. Прямой решатель имеет варианты для систем с распределённой памятью, при этом используется MPI. Кроме этого, существует версия для графических ускорителей, разработанная Acceleware Corporation. Позднее версия была портирована на гетерогенный кластер [8].
Прямой решатель может работать с симметричными и несимметричными матрицами, причём начиная с версии 6.13 графические ускорители могут быть использованы и для несимметричных матриц. Итерационный решатель использует метод подпространств Крылова с предобуславливателем [6]. Авторам не удалось найти информацию об ускорении итерационного решателя с помощью графических ускорителей.
5. DEFORM3D
Данный программный комплекс специально предназначен для моделирования больших пластических деформаций в инженерных задачах. Частичная информация об используемых в пакете решателях описана в руководстве пользователя [16].
КЭПК позволяет по отдельности выбирать метод для решения СЛАУ, возникающих в задачах теплопереноса и деформации. Для моделирования механики можно использовать прямой метод или методы CG и GMRES. Параллелизм достигается с помощью MPI.
6. QFORM3D
В руководстве пользователя [30], содержащем отдельный раздел по теоретическим основам решения задач предметной области, сообщается, что СЛАУ при моделировании деформаций решаются прямым методом с учётом разреженности матриц.
7. COMSOL Multiphysics
Несколько менее известным коммерческим программным комплексом конечно-элементного анализа является COMSOL Multiphysics. Данный конечно-элементный программный комплекс поставляется с несколькими решателями [35]. Поставляемые прямые решатели основаны на LU-разложении и взяты из нескольких известных библиотек: MUMPS, Intel PARDISO, а также из SPOOLES [36]. Все решатели могут работать на гомогенных системах как с общей, так и с распределённой памятью. Итерационные решатели основаны на распространённых итерационных методах решения СЛАУ, таких как CG, GMRES, BiCGStab и FGMRES.
8. Fidesys
Российская компания «Фидесис» была основана в 2009 году и начиная с 2012 года выпускает конечно-элементный программный комплекс CAE Fidesys. Данный комплекс содержит модуль для работы на современных кластерных суперкомпьютерах, включая гетерогенные [18]. В руководстве пользователя [4] сообщается о наличии двух решателей: прямого и итерационного. Прямой решатель основан на LU-разложении. Детали реализации решателей не раскрываются.
Следует отметить, что CAE Fidesys изначально проектировалась с учётом гетерогенных систем и может использовать по крайней мере один ускоритель на всех этапах конечно-элементного анализа.
9. CalculiX
CalculiX является бесплатным КЭПК с открытым исходным кодом. Информация об используемых решателях описана в руководстве пользователя [17]. КЭПК может использовать решатель SGI [34], оригинальную версию PARDISO, SPOOLES [36], TAUCS [39] и внутренний итерационный решатель.
10. Code_Aster
Code_Aster является бесплатным КЭПК с открытым исходным кодом. Команды КЭПК, по-видимому, являются транслитерацией французского языка. Документация к комплексу также написана на французском языке, однако существует машинный вариант перевода на английский [13]. Машинный перевод позволяет сделать выводы об используемых решателях. В частности, используются прямые решатели, включающие в себя разработанный авторами мультифронтальный решатель, основанный на LDLT-разложении, подпрограммы из библиотеки MUMPS и прямой решатель, основанный на формате хранения матрицы «SkyLine».
Набор итерационных решателей включает в себя предобусловленный метод сопряжённых градиентов из библиотеки MUMPS, а также решатели библиотеки PETSc [3]. В качестве предобуславливателя используется неполное LU-разложение с ограничением на ширину ленты (ILU(k)). Кроме этого, комплекс включает в себя решатель, называемый авторами гибридным.
11. Elmer
Конечно-элементный программный комплекс Elmer разрабатывается финским научным сообществом и является бесплатным программным обеспечением с открытым исходным кодом. Документации по решению СЛАУ посвящена отдельная глава руководства пользователя [31]. Используемый по умолчанию прямой метод основан на процедурах библиотеки LAPACK [9] для ленточных матриц. Также имеется интерфейс для решателя Umfpack, в данный момент входящего в набор SuiteSparse [37]. Данный решатель использует мультифронтальную LU-факторизацию. Кроме этого, КЭПК поддерживает библиотеки MUMPS, SuperLU и стандартную версию PARDISO.
Комплекс поддерживает обширный список итерационных методов, включая CG, GMRES, CGS, TFQMR, BiCGStab и BiCGStab(l), а также обобщённый метод сопряжённых невязок (англ. «Generalized Conjugate Residual», GCR). Предобуславливатели включают в себя предобуславливатель Якоби, предобуславливатели ILU(k), ILU(0) и ILUT.
Кроме этого пакет Elmer может использовать алгебраический и геометрический мультисеточный методы.
Заключение
Большинство КЭПК общего назначения имеют по меньшей мере две группы взаимодополняющих решателей. Первая группа включает прямые решатели, предназначенные для работы с относительно небольшими задачами, а также задачами, порождающими СЛАУ с высокими числами обусловленности. Вторая группа представляет собой набор итерационных решателей, основанных на подпространствах Крылова с различными предобуславливателями.
Решатели первой группы (прямые) основаны на разложении Холецкого, LDLT-разложении или LU-разложении в зависимости от характера порождаемой задачей матрицы. Реализации прямых методов (включая библиотечные) используют схожий набор решений, включая ограниченный список алгоритмов предобработки матрицы и использование внешней памяти для хранения факторизаций.
Итерационные решатели отличаются бо́льшим разнообразием. Распространены методы, основанные на подпространствах Крылова. Несмотря на то, что в научной литературе представлено значительное количество итерационных методов, подавляющее большинство производителей КЭПК используют методы, рекомендуемые пособием [14], при этом методы CG, BiCGStab и GMRES используются значительно чаще, чем метод CGS и метод минимальных невязок (англ. «Minimal Redidual», MinRES). Вероятно, это объясняется универсальностью используемых в КЭПК методов. Однако сходимость итерационных методов сильно зависит от параметров матрицы, при этом даже на последовательном вычислителе нельзя выделить однозначно подходящий для большинства случаев алгоритм.
Наиболее широко используются алгебраические предобуславливатели, в основном варианты неполного LU-разложения. Другие предобуславливатели используются значительно реже.
В КЭПК зачастую используются готовые бесплатные и коммерческие библиотеки подпрограмм для решения СЛАУ прямыми методами, описанные в четвёртом разделе данной статьи. Библиотеки, содержащие итерационные методы, такие как PETSc, включаются реже. Можно предположить, что авторы библиотек конструируют итерационные методы с помощью библиотек математических примитивов. Вероятно, со временем в поставку КЭПК будут включаться новые коммерческие итерационные решатели.
Ведущие разработчики программного обеспечения уделяют значительное внимание адаптации решателей под новые архитектуры высокопроизводительных систем, включая вычислительные ускорители. Наиболее часто применяемым стандартом программирования ускорителей является платформа Nvidia CUDA. Её используют КЭПК ANSYS, Abaqus FEA и FIDESYS. Комплекс ANSYS также поддерживает гетерогенные системы с Intel Xeon Phi.
Работа выполнена в рамках Программы № 7 УрО РАН, проект № 15-7-1-17.
Рецензенты:
Суханов В.И., д.т.н., доцент, заведующий кафедрой программных средств и систем факультета ускоренного обучения, ФГАОУ ВПО «Уральский федеральный университет имени первого Президента России Б.Н. Ельцина», г. Екатеринбург;
Залазинский А.Г., д.т.н., профессор, заведующий лабораторией системного моделирования, ФГБУН «Институт машиноведения» Уральского отделения Российской академии наук», г. Екатеринбург.