CC BY
52
Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2010, № 1, с. 202-209
УДК 519.612
ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ НА СКИФ МГУ «ЧЕБЫШЕВ» КОМБИНИРОВАННОЙ MPI+THREADS-РЕАЛИЗАЦИИ АЛГОРИТМА РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ, ВОЗНИКАЮЩИХ ВО FlowVision ПРИ МОДЕЛИРОВАНИИ ЗАДАЧ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ*
© 2010 г. Г.Б. Сушко, С.А. Харченко
ООО «ТЕСИС», Москва skh@tesis.com.ru
Поступила в редакцию 14.04.2009
Моделирование задач вычислительной гидродинамики с использованием неявных методов приводит к системам линейных уравнений, описываемым неструктурированными разреженными матрицами большой размерности. Во FlowVision при их решении используются предобусловленные итерационные алгоритмы типа подпространства Крылова с предобусловливателем высокого порядка точности. Для уменьшения времени вычислений этот алгоритм распараллеливается по распределенной памяти на основе технологии MPI, а по общей памяти на основе Intel® Threading Building Blocks. Для явного выделения параллелизма в зависимых вычислениях используется упорядочивание типа вложенных сечений Nested Dissection. Приводятся результаты численных экспериментов по масштабируемости предложенных алгоритмов на СКИФ МГУ “Чебышев”.
Ключевые слова: система линейных уравнений, итерационые методы, неполная треугольная факторизация, параллельные вычисления, MPI, общая память, потоки вычислений.
Введение
Моделирование задач вычислительной аэро-и гидродинамики является актуальной проблемой многих отраслей промышленности. Параллельная версия программного комплекса Р1о"^ У1в1оп, разрабатываемого в ООО «ТЕСИС», используется для моделирования широкого круга индустриальных задач, включая задачи с подвижными телами и свободными поверхностями для областей со сложной геометрией. Для адекватного воспроизведения тонких физических эффектов в геометрически сложных трехмерных областях требуются подробные расчетные сетки, содержащие от сотен тысяч до десятков и даже сотен миллионов расчетных ячеек. Подобные задачи требуют огромных вычислительных ресурсов и могут быть решены только на самой современной мощной параллельной вычислительной технике [1].
Использование неявных численных методов во F1owVision приводит к необходимости решения систем линейных уравнений
Ах = Ь С1)
* Статья рекомендована к печати программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2009» (http:// agora.guru.ru/pavt).
с неструктурированными сильно разреженными матрицами большой размерности. В настоящее время для эффективного решения таких систем линейных алгебраических уравнений (СЛАУ) используются предобусловленные итерационные методы. Во F1owVision в качестве предо-бусловливателя используется ШШ, являющийся несимметричным обобщением предо-бусловливателя 1СН2 [2]. В работе [2] доказаны существование и устойчивость этого неполного разложения для произвольных симметричных положительно-определенных матриц, а также получены оценки скорости сходимости метода сопряженных градиентов с этим предобусловливанием. На практике пре-добусловливатели 1СН2 и ШШ показали свою высокую надежность и эффективность при решении разнообразных систем линейных уравнений. В качестве итерационной схемы используются методы типа подпространства Крылова: метод сопряженных градиентов СО
[3] в симметричном случае, алгоритм ОМКБ8
[4] и алгоритм Ланцоша [3] в несимметричном случае.
Для уменьшения времени счета проводимые вычисления требуется распараллелить в соответствии с архитектурой используемой вычислительной техники. Доступная в настоящий момент техника обладает как правило архитек-
турой с неоднородным доступом к памяти. А именно, с одной стороны, имеется набор вычислительных узлов с распределенной памятью, обмен данными между которыми может быть осуществлен по быстрой обменной сетке. С другой стороны, каждый узел представляет собой многопроцессорный/многоядерный компьютер с общим доступом к оперативной памяти.
В предыдущих работах [5-8] авторами были предложены алгоритмы распараллеливания вычислений при решении систем уравнений для случая параллельных вычислений с распределенной [5-7] на основе MPI и с общей [8] памятью на основе OpenMP/TBB (ТВВ - Intel Threading Building Blocks) [9, 10]. Алгоритмы, описанные в этих работах, основаны на единых подходах к распараллеливанию вычислений, базируются на рекурсивной геометрической декомпозиции расчетной сетки на основе упорядочиваний типа вложенных сечений Nested Dissection (ND) [11], используют специальные подходы к распараллеливанию вычислений с поверхностными межпроцессорными границами [8] и при этом учитывают адаптацию расчетной сетки [7].
Эффективная хорошо масштабируемая организация параллельных вычислений при большом общем числе доступных процессоров/ядер предполагает хорошую балансировку вычислений на различных процессорах между собой. Чем больше число задействованных процессоров/ядер, тем более точная балансировка вычислений требуется для поддержания заданного уровня масштабируемости. К сожалению, в случае использования распределенной памяти принимать решение об организации вычислений (о распределении данных по процессорам и
о взаимодействии процессоров/ядер) приходится до собственно самих вычислений. Это связано с тем, что динамический перенос вычислений с одних процессоров на другие в случае распределенной памяти требует слишком больших обменов необходимыми данными. При этом в случае использования распределенной памяти трудно хорошо сбалансировать вычисления, если объем этих вычислений заранее точно неизвестен. С другой стороны, если вычисления проводятся по общей памяти, то имеется возможность динамически определить, какое ядро (процессор) должно проводить текущее вычисление, так как можно организовать вычисления по принципу клиент - сервер с потенциально наилучшей возможной балансировкой нагрузки.
В работе описывается параллельная реализация алгоритма решения СЛАУ с предобуслов-
ливателем ICH2/ILU2 и итерационной схемой типа подпространства Крылова для компьютеров с неоднородным доступом к памяти. Для уменьшения времени вычислений этот алгоритм заранее распараллеливается по распределенной памяти (по узлам вычислений) с использованием стандарта MPI, а по общей памяти узла распараллеливание по процессорам/ядрам осуществляется с динамическим распределением нагрузки на основе технологии Intel® Threading Building Blocks [10]. Приводятся результаты численных экспериментов по масштабируемости предложенных алгоритмов на СКИФ МГУ «Чебышев» для систем линейных уравнений, возникающих при моделировании во FlowVision задач вычислительной гидродинамики.
1. Блочные биения и матрица зависимостей блочных вычислений при решении систем уравнений
Естественным способом явного выделения параллелизма вычислений является введение для основных структур данных подходящей крупноблочной структуры. При этом явное указание зависимостей вычислений при построении объектов этой блочной структуры естественным образом приводит к понятию матрицы зависимостей блочных вычислений. В терминах такой матрицы для каждого объекта крупноблочной структуры можно указать, от каких других объектов крупноблочной структуры этот объект зависит. Тем самым в параллельных вычислениях возникает частичный порядок блочных вычислений.
В данном разделе вводится понятие блочного биения и блочной структуры матрицы системы линейных уравнений, а также в терминах блочных операций описываются основные алгебраические вычисления алгоритма решения системы уравнений. Вводится понятие матрицы зависимостей блочных вычислений, описывается связь этой матрицы со структурой параллельных вычислений для основных вычислений алгоритма. Также в этом разделе приводится способ построения упорядочивания и построения блочного биения для системы линейных уравнений, для которого имеет место большой ресурс параллелизма вычислений и для которого матрица зависимостей блочных вычислений обладает заданными свойствами.
1.1. Упорядочивание и блочное биение
Предположим, что все неизвестные системы уравнений пронумерованы некоторым образом.
Нумерация переменных может быть различной в различных случаях и для различных целей. Например, нумерация может быть специально подобрана так, чтобы в ней явным образом выделялся параллелизм вычислений для заданного числа процессоров и/или чтобы вычисления были максимально эффективны при решении системы уравнений, и т.д. Для выбранной нумерации неизвестных системы уравнений вводим разбиение множества всех индексов неизвестных на некоторое количество подмножеств, например р , р > 1, каждое из которых содержит набор индексов неизвестных, расположенных подряд. Нумеруем выбранные множества индексами от 1 до р в соответствии с возрастанием индексов в множествах. Пусть матрица системы уравнений А переупорядочена в соответствии с новой нумерацией неизвестных. Обозначим через Лу блок матрицы системы
уравнений, отвечающий наборам индексов I по строке и у по столбцу. Матрица системы уравнений в терминах введенного блочного биения имеет вид:
А =
Аіі - Аір
1 > а12 0 0 11 «12 1
а21 а22 а23 = ^21 ^22 0 0 и22 «23
1 а32 ^32 0 0 1
(3)
Здесь блок а22 имеет размер 1, блоки а21, а23, 121, и23 - строки, а а12, а32,132, и12 - столбцы. Рекуррентный обход матрицы производится последовательно по строкам V и столбцам Ь сверху вниз, т.е. блоки, Ьп, иі, і, у = 1, к ,3, уже
вычислены, блоки 122, /32, и22, и23 рассчитываются на данном шаге, а блоки £33 и изз пока неизвестны. На каждом шаге вычислений необходимо:
1. Из уравнения 122и22 = а22 - 121ип вычислить значения /22, и22.
(2)
Для перенумерованного вектора правой части Ь обозначим через Ъ{ блок этой правой части, соответствующий набору индексов I из (1). Для остальных векторных и матричных данных алгоритма решения системы уравнений (1) их блочные биения вводятся аналогичным образом.
1.2. Расширенная блочная разреженность как матрица зависимостей блочных вычислений
Алгоритм решения систем уравнений, используемый во FlowVision, предполагает неполную факторизацию матрицы системы уравнений и последующие предобусловленные итерации. Построим для этих вычислений описание максимально возможной зависимости по данным с учетом разреженности матрицы системы уравнений.
Неполная треугольная факторизация
Следуя изложению работы [8], рассмотрим алгоритм (точной) треугольной факторизации -разложения матрицы А в произведение нижней Ь и верхней и треугольных матриц. В блочном виде рекуррентное соотношение этой операции можно записать следующим образом:
2. Из уравнений и23 = (а23 - 121и13 )/22 и
132 = (а23 - Ь31и12 )и221 найти значения и23 и 132.
Различные варианты неполных разложений, в том числе и 1СИ2/1Ьи2 [2], получаются из описанного выше алгоритма путем некоторой модификации шагов 1 и 2.
Зависимость по данным в этом алгоритме для текущей строки и показана на рис. 1, для столбца Ь зависимость по данным аналогична. Таким образом, при вычислении текущей строки и/столбца Ь используются некоторые ненулевые значения из уже посчитанных строк и/столбцов Ь.
Рис. 1. Зависимость по данным при факторизации строки и
Неполная треугольная факторизация в терминах блочных вычислений. Матрица зависимостей блочных вычислений
В терминах блочного биения (2) неполные вычисления неполной факторизации могут быть снова записаны в виде шагов 1 и 2 с той лишь разницей, что количество шагов теперь будет равным числу блоков блочного биения матрицы и что в новых формулах вычислений операция умножения элементов матриц заменяется на операцию умножения разреженных подматриц, а операция умножения на обратное диагональное значение заменяется в зависимости от вари-
анта неполного разложения соответствующим приближенным блочным решением систем уравнений с диагональным блоком с учетом разреженности участвующих подматриц. Зависимости по данным при этих вычислениях в терминах блоков снова описываются диаграммой, показанной на рис. 1, при этом зависимости по данным имеют место, только если соответствующие блоки матриц не равны тождественно нулю.
Рассмотрим вопрос о блочной структуре разреженности матриц предобусловливания Ь и и. С точки зрения качества предобусловлива-ния, очевидно, блочное заполнение матриц Ь и и должно быть не меньшим, чем блочное заполнение исходной матрицы системы уравнений. С другой стороны, как это показано выше, максимальное блочное заполнение матриц Ь и и предобусловливания соответствует так называемой символьной факторизации блочной разреженности матрицы системы уравнений. На практике, с учетом неполноты вычислений блочная разреженность матриц Ь и и занимает некоторое промежуточное положение между этими двумя естественными пределами. Алгоритм распараллеливания вычислений с поверхностными межпроцессорными границами, описанный в работе [7], также предполагает некоторое дополнительное уменьшение заполнения блочной разреженности матриц Ь и и по сравнению со случаем символьной факторизации блочной структуры разреженности матрицы системы уравнений.
С учетом вышесказанного, для обеспечения фиксированной блочной зависимости параллельных вычислений и для обеспечения гарантированного достаточного ресурса параллелизма при вычислениях для каждого блочного биения матрицы будем рассматривать матрицу зависимостей блочных вычислений. Это матрица, структура разреженности которой есть сумма допустимых с учетом всех необходимых требований блочных разреженностей матриц Ь и и предобусловливания. По определению будем предполагать, что структура разреженности матрицы зависимостей блочных вычислений не меньше, чем блочная структура разреженности матрицы системы уравнений, и не больше, чем сумма разреженностей факторов символьной факторизации матрицы блочной разреженности системы уравнений.
В терминах матрицы зависимостей блочных вычислений легко описываются все зависимости по данным, которые имеют место при построении неполных факторов Ь и и. С другой стороны, соответственно нижняя и верхняя тре-
угольные части этой же матрицы могут быть использованы для описания зависимостей блочных вычислений при решении системы уравнений с нижней Ь и верхней и треугольной матрицей. В силу вложения блочной структуры разреженности матрицы системы в структуру разреженности матрицы зависимостей блочных вычислений эта же матрица может быть использована при описании зависимостей блочных вычислений при умножении матрицы системы уравнений на вектор. Получается, что структура разреженности матрицы зависимостей блочных вычислений тем или иным образом описывает зависимости по данным на всех этапах алгебраического алгоритма решения системы линейных уравнений и вводит частичный порядок блочных вычислений.
1.3. Нумерация ячеек, блочные биения и декомпозиция для параллельных вычислений
Эффективная организация параллельных вычислений требует выделения независимых вычислений, которые можно было бы производить одновременно, а также зависимых вычислений, которые можно было бы выполнять максимально сбалансированно и с минимальными задержками. Для того чтобы выделить подобные вычисления, при решении систем уравнений во FlowVision в случае использования адаптированных расчетных сеток используются подходы, описанные в работах [5-7].
Сначала производится объемная декомпозиция расчетной сетки в терминах так называемых начальных ячеек [6]. Каждая начальная ячейка представляет собой элемент регулярной декартовой начальной сетки, который может содержать некоторое множество реальных ячеек. Количество реальных ячеек в начальных ячейках увеличивается при адаптациях расчетной сетки в соответствии с заданными пользователем критериями. Объемная декомпозиция осуществляется на число процессоров/ядер -степень двойки - и производится последовательным биением на две подобласти каждого предыдущего разбиения. Объемная декомпозиция на две подобласти осуществляется с помощью программы РагМе^8 [12] с учетом условий минимального возможного числа реальных ячеек на границе между подобластями и с учетом условия наилучшего возможного совпадения полного числа расчетных ячеек в обеих подобластях.
На основе объемной декомпозиции задачи строятся однослойные поверхностные границы между подобластями [7], эти границы по построению проходят по поверхностным грани-
цам начальных ячеек. Внутри каждой полученной объемной подобласти ячейки нумеруются в соответствии с упорядочиванием вложенных сечений ND [11]. Для ячеек на границе между подобластями при помощи специального алгоритма предсказывается разреженность дополнения по Шуру для этой границы [7]. На основе предсказанной разреженности дополнения по Шуру производится декомпозиция поверхностных границ на близкие по числу ячеек поверхностные подобласти по принципу непосредственного примыкания этих поверхностных подобластей к исходным объемным подобластям. Это сделано для выделения максимального параллелизма при работе с поверхностными границами и для минимизации обменов и зависимостей в вычислениях при их обработке. Возникающая при каждой поверхностной декомпозиции одномерная граница между поверхностными подобластями ассоциируется с одной (любой из примыкающей пары) объемной подобластью. После перенумерации матрицы в соответствии с таким образом построенным объемным и поверхностным биением возникают структуры разреженности, показанные для тестовой задачи на рис. 2. На рис. 3 для той же тестовой задачи показаны соответствующие матрицы зависимостей блочных вычислений.
2. MPI+Threads-реализация алгоритма решения систем уравнений на компьютерах с неоднородным доступом к памяти
На компьютерах с неоднородным доступом к памяти возможны различные способы организации вычислений. В данной работе при распараллеливании вычислений между узлами для организации обменов и синхронизаций используется стандарт обмена сообщениями MPI. В каждом узле с общей памятью при распараллеливании по процессорам/ядрам для организации синхронизаций по общей памяти используется технология Intel® Threading Building Blocks [10]. При этом предполагается, что на каждом узле имеется только один MPI-процесс, который затем порождает на этом узле нужное количество одновременно работающих потоков вычислений.
При работе с распределенной и с общей памятью естественными являются различные подходы к распараллеливанию вычислений. С одной стороны, на распределенной памяти для минимизации обменов естественно заранее распределить данные по процессорам и при этом по возможности сбалансировать вычисления. С другой стороны, при работе по общей памяти
нет необходимости заранее распределять данные по процессорам/ядрам. Для наилучшей возможной балансировки вычислений можно динамически выбирать ядро, которое будет проводить текущие вычисления.
На компьютерах с неоднородным доступом к памяти для достижения по возможности наилучших результатов естественно использовать одновременно оба подхода к распараллеливанию вычислений. А именно, сначала все вычисления заранее распределяются по узлам таким же образом, как это делалось при работе только с распределенной памятью. Затем каждое вычисление в узле динамически распараллеливается по процессорам/ядрам с общей памятью таким же образом, как это делается в случае только общей памяти. При этом для достижения эффекта самобалансировки вычислений внутри узла полное число объектов, с которыми предполагается проводить вычисления внутри одного узла, должно быть достаточно большим. Это необходимо для обеспечения возможности организации вычислений по принципу клиент/сервер с динамическим выбором ядра, на котором будет проводиться текущее вычисление.
Описанная методика организации параллельных вычислений на компьютерах с неоднородным доступом к памяти в применении к алгоритмам решения систем уравнений реализована следующим образом. Пусть имеется з узлов и ? процессоров/ядер в каждом узле, итого всего р = 5? процессоров/ядер. Производится параллельная декомпозиция задачи (в том числе построение блочного биения матрицы системы уравнений и матрицы зависимостей параллельных вычислений) на число процессоров/ядер, равное pq, где д > 1 - целое число, обозначающее запас параллелизма для автоматической балансировки вычислений внутри узла. В данной работе использовалось значение д = 4 . Далее, на каждый узел распределяются все данные, соответствующие tq условным процессорам/ядрам этого узла. Все параллельные вычисления проводятся в соответствии с блочными структурами соответствующих матриц на основе матрицы зависимостей блочных вычислений, как это описано в разделе 1.2.
3. Результаты численных экспериментов
Для проведения численных экспериментов была выбрана задача об обтекании спортивного автомобиля. Адаптивная расчетная сетка для этой задачи содержит 2 994 497 расчетных
Рис. 2. Разреженность матрицы системы уравнений с учетом блочного биения для 16 и 64 процессоров/ядер
Рис. 3. Матрица зависимостей блочных вычислений для 16 и 64 процессоров/ядер
Таблица
Результаты численных экспериментов
Конфигурация 2 узла по 4 ядра 4 узла по 2 ядра 4 узла по 4 ядра 4 узла по 8 ядер
Алгоритм МРІ МРІ+ТББ МРІ МРІ+ТББ МРІ МРІ+ТББ МРІ+ТББ
Факторизация 14.4 9.32 9.71 9.16 7.63 5.01 3.98
Итерации 102.27 54.13 51.69 54.39 53.08 29.14 31.74
ячеек. Структура адаптаций сетки в различных ракурсах показана на рисунке 4.
Объемная декомпозиция расчетной сетки на 16 процессоров в терминах начальных ячеек [6] показана на рисунке 5. На рисунке 2 показана разреженность матрицы коэффициентов системы с учетом блочного биения для случаев 16 и 64 процессоров/ядер, построенная в соответствии с описанием раздела 1.3. На рисунке 3 представлены структуры разреженности матрицы зависимостей блочных вычислений для случаев 16 и 64 процессоров/ядер.
Численные эксперименты проводились на СКИФ МГУ «Чебышев». В таблице представлены результаты численных экспериментов.
Из результатов численных экспериментов следует, что комбинированная MPI+Threads-версия алгоритма решения систем уравнений позволяет ускорить вычисления до 2 раз по сравнению с чистой МР1-версией. Это связано с сильным дисбалансом вычислений, который естественным образом возникает при использовании неполных факторизаций по значениям в условиях существенной адаптации расчетной
сетки. Этот дисбаланс удается компенсировать за счет динамической балансировки вычислений внутри узла. При этом положительный эффект наблюдается начиная с 4 ядер в узле, при использовании 2 ядер на узел в итерационной схеме наблюдается даже небольшой проигрыш по времени по сравнению с чистой МР1-версией. Из таблицы также видно, что увеличение числа используемых ядер на узел до 8 привело к дальнейшему ускорению неполной факторизации и к небольшому замедлению итерационной схемы. Этот эффект, по-видимому, связан с особенностями переиспользования памяти в этих алгоритмах. В неполной факторизации каждое данное переиспользуется множество раз, в отличии от итерационной схемы, а потому неполная факторизация слабее подвержена влиянию эффекта недостаточной масштабируемости узла с точки работы с оперативной памятью при большом числе используемых ядер.
Заключение
В работе представлен комбинированный MPI+Threads параллельный алгоритм решения системы линейных уравнений для компьютеров с неоднородным доступом к памяти. На примере системы линейных уравнений, возникающей во FlowVision при моделировании тестовой задачи вычислительной гидродинамики, показана его эффективность по сравнению с чистой MPI-реализацией.
Список литературы
1. Аксенов А.А., Дядькин А.А., Кутин В.А. и др. Решение больших задач вычислительной гидродинамики на СКИФ МГУ с помощью FlowVision // Материалы Всероссийской научной конференции "Научный сервис в сети Интернет: решение больших задач", Новороссийск, 22-27 сентября 2008 г. Изд-во Московского университета, 2008. С. 69-73.
2. Kaporin I.E. High Quality Preconditioning of a General Symmetric Positive Definite Matrix Based on its
UTU + UтR + RTU decomposition // Numer. Linear Algebra Appl. 1998. V. 5. P. 483-509.
3. Тыртышников Е.Е. Краткий курс численного анализа. М.: ВИНИТИ, 1994.
4. Saad Y., Schultz M.H. GMRES: A generalized minimum residual algorithm for solving non-symmetric linear systems // SIAM J. Sci. Comput. 1986. 7. Р. 856-869.
5. Харченко С.А. Параллельная реализация алгоритма решения систем линейных уравнений в пакете FlowVision // Труды конференции «Инженерные системы-2007». 2007. С. 135-144.
6. Дядькин А.А., Харченко С.А. Алгоритмы декомпозиции области и нумерации ячеек с учетом локальных адаптаций расчетной сетки при параллельном решении систем уравнений в пакете FlowVision // Труды Всероссийской научной конференции «Научный сервис в сети Internet: многоядерный компьютерный мир». 2007. С. 201-206.
7. Харченко С.А. Влияние распараллеливания вычислений с поверхностными межпроцессорными границами на масштабируемость параллельного итерационного алгоритма решения систем линейных
уравнений на примере уравнений вычислительной гидродинамики // Труды Международной научной конференции «Параллельные вычислительные технологии» (ПаВТ'2008). Санкт-Петербург, 28 января -
1 февраля 2008 г. Челябинск: Изд. ЮУрГУ, 2008. С. 494-499.
8. Сушко Г.Б., Харченко С.А. Многопоточная параллельная реализация итерационного алгоритма решения систем линейных уравнений с динамическим распределением нагрузки по нитям вычислений // Труды Международной научной конференции «Параллельные вычислительные технологии» (ПаВТ'2008). Санкт-Петербург, 28 января - 1 февраля 2008 г. Челябинск: Изд. ЮУрГУ, 2008. С. 452457.
9. OpenMP Application Program Interface - 2.5 // OpenMP Architecture Review Board. 2005.
10. Intel Threading Building Blocks Tutorial - 1.6 // Intel Corp., 2007.
11. George A., Liu J.W. Computer Solution of Large Sparse Positive Definite Systems. Prentice Hall, 1981.
12. http://glaros.dtc.umn.edu/gkhome/metis/parmetis/ overview
EXPERIMENTAL STUDY ON SKIF MSU “CHEBYSHEV” OF THE COMBINED MPI+THREADS IMPLEMENTATION OF THE ALGORITHM FOR SOLVING SYSTEMS OF LINEAR EQUATIONS ARISING IN FLOWVISION WHEN SIMULATING COMPUTATIONAL FLUID DYNAMICS PROBLEMS
G.B. Sushko, S.A. Kharchenko
Simulation of CFD problems in FlowVision by implicit numerical methods requires the solution of systems of linear equations described by large unstructured sparse matrices. In FlowVision the solution algorithm for such systems of linear equations includes the computation of the high order incomplete triangular factors (ICH2/ILU2) of the coefficient matrix followed by calculations using the algorithm of the preconditioned Krylov subspace iterative method. In order to reduce the solution time, the algorithm is a priori parallelized over the distributed memory using the MPI standard, while over the common memory it is parallelized by using Intel® Threading Building Blocks (TBB) with automatic dynamic load balancing. For extracting explicit parallelism in dependent computations and for solving parallel systems of equations with triangular matrices, the nested dissection (ND) ordering algorithm is used. The numerical experiment results on scalability of the proposed algorithms performed on SKIF MSU “Chebyshev” demonstrate their efficiency in simulating CFD problems in FlowVision. The results of numerical experiments performed on SKIF MGU “Chebyshev” demonstrate the efficiency and parallel scalability of the proposed algorithms when simulating sample CFD problems in FlowVision.
Keywords: system of linear equations, iterative methods, incomplete triangular factorization, parallel computations, MPI, common memory, computational threads.
