Многочисленные исследования вычислительных схем на базе метода конечных элементов и метода граничных элементов привели к выводу про их формальное родство со многими аппроксимационными методами. Если рассматривать все эти методы из позиций вариационных принципов, то можно исследовать все варианты реализации вычислительного процесса. Использование свойств математического языка, описание вариационных функционалов, вариантов дальнейших вычислений создало возможность управлять вычислительным процессом с помощью формульного языка, на которой описываются простые формулы и выводы функционалов. Формульный язык разрешает в рамках единого вариационного формализма решать и исследовать задачи математической физики в разной вариационной постановке и с разным функциональным базисом.
Общая вариационная постановка задачи математической физики предполагает наличие некоторой области, некоторого списка базисных величин ui, которые варьируются в заданной области, и некоторого функционала F(u1,u2...,x,t) для области, что необходимо минимизировать. Дополнительно необходимое наличие предельных и, если это динамическая задача, начальных условий.
Большинство существующих численных методов решения краевых задач основано на аппроксимации решения с помощью определенной системы координатных (базисных) функций: Ui = сумма по Cij, причем предполагается полнота системы координатных функций в пространстве допустимых решений и сходимость линейной комбинации координатных функций к точному решению по мере увеличения числа ее членов. Сюда можно отнести такие широко используемые средства вычислительной математики, как метод конечных элементов, метод предельных элементов, метод Ритца, метод наименьших квадратов, метод коллокаций, расчет балок и пластин с помощью собственных функций соответствующего дифференциального оператора. Общая задача определения минимума функционала сводится при этом к решению системы алгебраических уравнений относительно неизвестных коэффициентов Ci.
Однако подобная методика накладывает и некоторые важные ограничения на вид функционалов, которые минимизируются. Любой функционал, минимизируемый по данной методике, должен соответствовать краевой задаче с линейным дифференциальным оператором. Поэтому все многообразие функционалов сводится к поверхностным и объемным интегралам от квадратичной формы варьируемых величин. Впрочем, в квадратичной форме могут принимать участие как коэффициенты весовые функции, но они должны быть известны или быть вычислены раньше.
Чтобы войти в подынтегральное выражение, любая формула для величины, которая варьируется, должна представлять собой линейную комбинацию координатных функций. Поэтому формулы для варьируемых величин должны строится с помощью других варьируемых величин и базисных варьируемых величин только с помощью операций сложения, вычитания, умножения деления на константу, дифференцирования.
Разработан стандарт языка FORTU, предназначенного для описания схем решения краевых задач математической физики. Существующий стандарт языка есть обобщением результатов продолжительного развития инструментальных систем, разрабатываемых авторским коллективом во главе с профессором Толоком В.А.
Возникновение МКЭ связано с невозможностью аналитического решения большинства задач механики деформированного твердого тела и механики разрушения. При решении задач МКЭ исследуемый объект представляется в виде совокупности дискретных подобластей - конечных элементов. Внутри каждого конечного элемента искомая непрерывная величина аппроксимируется совокупностью кусочно-непрерывных функций. Искомая величина внутри КЭ определяется с помощью значений этой величины в конечном числе точек исследуемого объекта - узловых точках, или узлах, которые, как правило, являются характерными точками КЭ.
При решении задачи искомая величина подчиняется совокупности аналитических зависимостей, которые моделируют свойства материала (непрерывность, изотропность и др.) и законы механики (закон равновесия сил, вариационные принципы). Это позволяет при заданных начальных и граничных условиях получить однозначное решение задачи.
Система FORTU - новый инструмент, пригодный для решения современных задач конструирования и компьютерного моделирования. Существующие программные комплексы для расчета задач механики деформируемого твердого тела - COSMOS, ANSYS, и др. - нацелены на массовое решение типичных задач с определенной математической моделью. Переход к другой математической модели есть для них часто невозможным. Но конструирование само по себе не является стандартизированным тривиальным процессом.
Система FORTU направлена именно на то, чтобы дать инженеру возможность уточнять и видоизменять математическую модель в процессе конструирования. Она может послужить дополнением к существующим программным комплексам на поисковом этапе конструирования. Система FORTU реализует транслятор одноименного языка для построения векторов и разреженных матриц и описания схемы манипуляции этими объектами для получения численного решения.
При разложении разреженной матрицы A общего вида для обеспечения численной устойчивости необходима та или иная форма выбора главного элемента, т.е. перестановки строк и/или столбцов. Таким образом, при заданной матрице A обычно получают разложение для P*A*QT, где P и Q - матрицы перестановок соответствующих размеров. Заметим, что умножение на P слева переставляет строки A, а умножение на Q справа переставляет столбцы A.
При разложении разреженной матрицы общего вида эти перестановки определяются в процессе разложения, путем компромисса между конкурирующими требованиями численной устойчивости и разреженности. Различные матрицы, хотя бы они и имели одинаковую структуру нулей-ненулей, обычно приводят к различным P и Q и, следовательно, имеют множители с различной структурой разреженности. Другими словами, для разреженных матриц общего вида, как правило, нельзя предсказать, где произойдет заполнение, пока не начались собственно вычисления. Это вынуждает пользоваться различными схемами динамического хранения, в которых память для заполнения выделяется в ходе разложения.
Матрицы, получаемые в инженерно-конструкторских и научных вычислениях, обычно соответствуют минимуму некоторого энергетического функционала и являются симметричными, положительно определенными. Симметричное гауссово исключение в применении к симметричной положительно определенной матрице не требует перестановок для поддержания численной устойчивости. Прямые методы для таких матриц основаны на разложении матрицы A в произведение треугольных матриц и требуют хранения треугольного множителя (целиком или частично). Метод Холецкого A=L*LT, где L - нижняя треугольная матрица, применим для положительно определенных матриц. Метод Мартина A=L*D*LT, где D - дополнительная диагональная матрица, применим как для положительно определенных, так и для более общих знаконеопределенных несингулярных эрмитовых матриц.
Поскольку P*A*PT также симметрична и положительно определена при любой матрице перестановки P, это значит, что можно симметрично переупорядочить, во-первых, не заботясь о численной устойчивости и, во-вторых, до начала реального численного разложения.
Эти возможности, обычно отсутствующие в случае матрицы A общего вида, имеют важнейшие практические последствия. Раз упорядочение можно определить до начала разложения, то можно определить также и местоположение заполнения, которое произойдет при разложении. Поэтому способ хранения L можно выбрать до реального численного разложения, так же как и зарезервировать место для элементов заполнения. Вычисления затем происходят при структуре хранения, остающейся статичной (неизменной).
Таким образом, три задачи: 1) выбор надлежащего упорядочения; 2) формирование подходящей схемы хранения; 3) реальные вычисления - могут быть разделены как самостоятельные объекты исследования и как разные модули программного обеспечения. Эта независимость задач имеет ряд отчетливых преимуществ. Она поощряет модульность при составлении математического обеспечения. Во многих инженерно-конструкторских приложениях приходится решать большое количество различных задач с положительно определенными матрицами, имеющими одинаковую структуру. Ясно, что упорядочение и формирование схемы хранения можно выполнить лишь однажды; поэтому желательно, чтобы эти этапы были изолированы от реальных вычислений.
Трудно получить аналитические выражения для числа арифметических действий при разложении и для размера полученного треугольного множителя (в терминах числа вхождений, отличных от нуля) для общей разреженной матрицы. Это потому, что вычисление и заполнение позиций ненулями в ходе разложения разреженной матрицы - функция от числа и позиций ненулей в первоначальной матрице.
Пусть A - заданная симметричная матрица, а P - матрица перестановки. Хотя структуры ненулевых элементов A и P*A*PT разные, их размеры одинаковы: Nonz(A) = Nonz(P*A*PT). Однако, и это обстоятельство большого значения, между Nonz(L(A)) и Nonz(L(P*A*PT)) для некоторой перестановки P разница может быть очень велика. В идеале мы хотели бы найти перестановку P*, которая бы минимизировала размер структуры ненулевых элементов матрицы заполнения: Nonz(L(P*A*PT)) = min Nonz(L(P*A*PT)).
До сих пор нет эффективного алгоритма для отыскания такой оптимальной P* в случае произвольной симметричной матрицы. Эта задача имеет сложность NP. Поэтому приходится прибегать к использованию эвристических алгоритмов, которые бы давали упорядочение с приемлемо малым, хотя и не обязательно минимальным размером заполнения Nonz(L(P*A*PT)).
В основе алгоритма - следующее наблюдение. Пусть узлы {x1, x2, ..., xi-1} исключены. Число ненулевых элементов в этих столбцах в дальнейшем не изменяется. Ясно, что для уменьшения числа ненулевых элементов i-го столбца в еще не разложенной матрице, на место i-го столбца нужно перевести столбец с наименьшим числом ненулевых элементов. Другими словами, схему можно рассматривать как метод уменьшения заполнение матрицы путем локальной минимизации заполнения для раскладываемой матрицы.
Процесс символического разложения матрицы можно описать в терминах теории графов как процесс рекурсивного построения последовательности графов исключения. Для получения очередного графа исключения необходимо удалить выбранный определенным образом узел, а потом превратить множество узлов, смежных исключаемому узлу, в заполненный подграф - клику. Граф заполнения состоит из ребер начального графа матрицы и ребер, добавленных при разложении. Мы должны знать этот граф, если намереваемся применить схему хранения, которая использует все нули треугольного множителя.
Разделение графа - важная проблема, которая имеет обширные применения во многих областях, включая научные вычисления, метод конечных элементов, линейное программирование, проектирование СБИС, транспортные задачи, планирование задач. Проблема состоит в том, чтобы разделить узлы графа на p приблизительно равных частей так, чтобы число ребер, соединяющих узлы из разных частей, было минимально.
Проблема разделения графа на k частей определяется следующим образом: имеем граф G=(V,E), |V|=n, делим V на k подмножеств, V1,V2,...,Vk так, что Vi AND Vj = 0 для i j, |Vi| приблизительно равно n/k, и сумма |Vi|=|V|, а число ребер из E, которые инцидентны узлам из различных подмножеств, минимально.
Проблема разделения графа на k частей может быть естественно расширена на взвешенные графы, которые имеют веса, ассоциированные с узлами и ребрами графа. В этом случае цель состоит в том, чтобы разделить узлы на k непересекающихся подмножеств так, что суммы весов узлов в каждом подмножестве приблизительно равны, а сумма весов ребер, которые инцидентны узлам из разных подмножеств, минимальна.
Разделение V на k частей обычно представляется вектором разделения P длины n, таким, что для каждого узла v V, P[v] - целое число между 1 и k, указывающее подмножество, которому принадлежит узел v. Для разделителя P число (или сумма весов) ребер, которые инцидентны узлам различных подмножеств, называется весом разделителя.
Граф G может быть разделен пополам, используя многоуровневый алгоритм. Основная структура многоуровневого алгоритма очень проста. Граф G сначала огрубляется до нескольких сотен узлов, вычисляется деление пополам этого намного меньшего графа, и затем это разделение проектируется (уточняется) назад к первоначальному графу. На каждом шаге восстановления грубого графа разделение улучшается (рафинируется). Так как более точный (больший) граф имеет большее число степеней свободы, такие обработки обычно уменьшают вес сепаратора.
Формально, многоуровневый алгоритм деления графа пополам работает следующим образом: рассмотрим взвешенный граф G0=(V0,E0), со взвешенными узлами и ребрами. Многоуровневый алгоритм деления графа пополам состоит из следующих трех стадий.
Алгоритмы вложенных сечений для минимизации заполнения при разложении разреженных матриц. Рекурсивное деление пополам применяется для поиска переупорядочения, сокращающего заполнение при разложении разреженной матрицы. Эти алгоритмы обычно упоминаются как алгоритмы вложенных сечений. Вложенные сечения рекурсивно разделяют граф на почти равные половины, удаляя узлы сепаратора, пока не получено желательное число разделений.
Один путь получения узлового сепаратора состоит в том, чтобы сначала получить деление пополам графа и затем вычислять узловой сепаратор из реберного сепаратора. Узлы графа нумеруются так, что на каждом уровне рекурсии, узлы сепаратора нумеруются после узлов в половинках графа. Эффективность и сложность алгоритма вложенных сечений зависит от алгоритма вычисления сепаратора. Вообще, маленькие сепараторы приводят к меньшему заполнению.
Большинство современных алгоритмов дискретизации трехмерных тел на КЭ базируются на идее средней оси. Средней ось (поверхностью) объекта называется геометрическое место точек являющихся центрами всех вписанных внутрь объекта окружностей (или сфер в трехмерном случае) максимального радиуса. Вписанная окружность максимального радиуса, по меньшей мере, должна касаться двух точек объекта. В противном случае - должна существовать другая вписанная окружность, полностью удовлетворяющая этим требованиям.
Современные генераторы КЭ-сети используют среднюю ось объекта для того, чтобы сформировать его подразбиение на 3-х - 6-ти сторонние односвязные области, каждая из которых в дальнейшем рассматривается как примитив, подлежащий дискретизации на КЭ.
Графическая библиотека OpenGL стала в последнее время промышленным стандартом при создании программ компьютерного моделирования и дизайна. Это гарантирует переносимость программного обеспечения и возможность применения аппаратных средств современных видеопроцессоров. При анализе полученного численного решения, как правило, возникает проблема наличия большого массива численной информации. Встроенная подсистема визуализации позволяет пользователю преодолеть эти трудности. Подсистема визуализации разрешает получить реалистическое изображение для исследуемой области с помощью цветовой градации, линий уровня и пр.
Для визуализации трехмерных объектов в отличие от существующих методик, которые акцентируют внимание на анализе поверхности, применяется методика сечений, которая дает возможность анализировать действительно "объемные" характеристики распределения. Сечение определяется как линия уровня f(x,y,z)=0, для которой определяется дискретное представление по конечно-элементному разбиению тела.
MPP (массивно-параллельная система) - система, состоящая из однородных вычислительных узлов, включающих: один или несколько центральных процессоров, локальную память (прямой доступ к памяти других узлов невозможен), коммуникационный процессор или сетевой адаптер, иногда - жесткие диски (как в SP) и/или другие устройства ввода/вывода. К системе могут быть добавлены специальные узлы ввода-вывода и управляющие узлы.
MPI - это программный инструментарий для обеспечения связи между ветвями параллельного приложения. MPI расшифровывается как "Message passing interface" ("Взаимодействие через передачу сообщений"). MPI предоставляет программисту единый механизм взаимодействия ветвей внутри параллельного приложения. MPI не зависит от машинной архитектуры: однопроцессорные или многопроцессорные системы с общей или раздельной памятью. MPI не зависит от взаимного расположения ветвей - на одном процессоре или на разных. MPI не зависит от API операционной системы (API = "applications programmers interface" = "интерфейс разработчика приложений"). Весь параллелизм явен: программист ответствен за правильное выделение параллелизма, и осуществление алгоритма, используя конструкции MPI.
Эффективное выполнение параллельных алгоритмов требует решения проблемы разделения графа, где узлы представляют вычислительные задачи, а ребра представляют межзадачный обмен данными. Узлам назначаются веса, пропорциональные количеству вычислений, выполненного каждой задачей, а ребрам - веса, пропорциональные количеству данных, которые должны быть переданы между задачами. Разделение этого графа вычислений на k частей может использоваться, чтобы распределить задачи между k процессорами. Так как разделение назначает на каждый процессор задачи, у которых суммарные веса одинаковы, то работа будет сбалансирована между k процессорами. Кроме того, поскольку алгоритм минимизирует вес сепаратора (ориентируясь на требование сбалансированной нагрузки), общий объем межзадачных коммуникаций также будет минимизирован.
Для разложения столбца, соответствующего какому-то неизвестному k, необходимо ранее разложить все столбцы, аккумуляция которых необходима для вычисления k-го столбца. Это соответствие порядка можно представить в виде древовидного графа, в котором отношения потомок-предок строятся на основе порядка разложения столбцов.
Разложение разреженной матрицы можно рассматривать как просмотр этого дерева исключения, при котором сначала раскладываются столбцы, соответствующие всем потомкам данного узла k, а затем строится разложение для самого данного узла k. В общем случае, для произвольной, в том числе и оптимальной перестановки неизвестных, дерево исключения не обладает какими-либо специальными свойствами. Оно может быть весьма произвольным и достаточно несбалансированным.
Анализ дерева исключения позволяет выявить в процессе разложения разреженной матрицы значительное количество вычислений, выполнение которых может быть произведено параллельно. Действительно, поддеревья одного корня дерева исключения могут быть вычислены одновременно. При параллельном разложении разреженной матрицы на первое место выходит не просто уменьшение количества требуемой памяти и объема необходимых вычислений, но создание такого дерева исключения, которое хорошо сбалансировано и имеет достаточно простой вид. Сбалансированность дерева исключения дает сбалансированность распределения вычислительной нагрузки на процессоры параллельной системы, а простота дерева - уменьшает сложность алгоритма назначения заданий процессорам.
В основе реализации системы памяти современных компьютеров лежат два принципа: принцип локальности обращений и соотношение стоимость/производительность. Принцип локальности обращений говорит о том, что большинство программ к счастью не выполняют обращений ко всем своим командам и данным равновероятно, а оказывают предпочтение некоторой части своего адресного пространства. Иерархия памяти строится из нескольких уровней, причем более высокий уровень меньше по объему, быстрее и имеет большую стоимость в пересчете на байт, чем более низкий уровень.
Адаптация методов матричной алгебры для таких систем предполагает переход от векторной формулировки к блочной формулировке. Наличие длинных вычислительных конвейеров требует специальных методов оптимизации для уменьшения числа зависимостей по данным и управлению в программах. Это достигается путем развертывания циклов.
МКЭ, FEM, FORTU, QMD, MMD, AMD, multilevel graph bisection, multilevel nested dissection, mesh generation, OpenGL, multifrontal sparse matrix factorization, MPI, MPP, BLAS.