CC BY
41
Известия Института математики и информатики УдГУ
2015. Вып. 2 (46)
УДК 519.63, 004.272.2
© А. К. Новиков, И. М. Кузьмин, Н. С. Недожогин
НЕКОТОРЫЕ ОСОБЕННОСТИ РЕАЛИЗАЦИИ КОНЕЧНО-ЭЛЕМЕНТНЫХ ТЕХНОЛОГИЙ ДЛЯ ГИБРИДНОЙ ВЫЧИСЛИТЕЛЬНОЙ АРХИТЕКТУРЫ1
Рассматриваются особенности параллельных конечно-элементных вычислений на гибридных архитектурах, связанные с заданием граничных условий первого рода (Дирихле). Сравниваются варианты задания граничных условий Дирихле при формировании и решении конечно-элементной системы уравнений, обеспечивающие сохранение симметрии матрицы коэффициентов.
Ключевые слова: параллельные вычисления, метод конечных элементов, граничные условия Дирихле, гибридные вычислительные архитектуры.
Введение
При конечно-элементном моделировании на кластерах гибридной (СРи+ОР11) архитектуры вычисления тем или иным образом распределяются между центральными и графическими процессорами. Соответственно, конечно-элементные данные, в том числе матрицы и вектора, располагаются на различных уровнях подсистемы памяти кластера: кэш-память, оперативная память, память ускорителя, подсистемы памяти других вычислительных узлов. Подобное размещение данных вносит особенности в алгоритмы, применяемые при реализации вычислительных схем метода конечных элементов (МКЭ) на гибридных архитектурах.
В методе конечных элементов применяется ряд алгоритмов, обеспечивающих сохранение симметрии матрицы коэффициентов системы сеточных уравнений при задании граничных условий первого рода или Дирихле. Такие алгоритмы применяются как при формировании конечно-элементных систем [1-3], так и при их решении [4].
Отметим, что изменения, связанные с граничными условиями Дирихле, возможно внести на нескольких уровнях: в собранной системе линейных алгебраических уравнений [1]; в локальных (элементных) матрицах жесткости и векторах [2], а также без изменения матрицы коэффициентов и вектора правой части [3].
Задание условий первого рода в параллельных конечно-элементных вычислениях на гибридной архитектуре требует модификации существующих алгоритмов. Далее рассматриваются особенности реализации граничных условий Дирихле на этапах формирования и решения конечно-элементных систем при параллельных вычислениях на гибридной архитектуре.
§ 1. Условия Дирихле на этапе формировании конечно-элементной СЛАУ
Рассмотрим, каким образом на гибридных архитектурах задаются граничные условия первого рода в случае явной сборки конечно-элементной системы уравнений. До задания граничных условий Дирихле имеем систему (1.1) с вырожденной матрицей коэффициентов К:
Ки = /; К = Кт; (Кх,х) = 0Ух = 0; К € и € / € (1.1)
Процесс сборки системы (1.1) запишем в следующем виде:
т т
К = £ СеТКеСе, / = £ СТ/е, (1.2)
е=1 е=1
здесь Се € ZWeXNV — матрица связи ости, Ке € — локальная матрица жесткости,
/е € — локальный вектор, N — число степеней свободы в конечном элементе е; т — число конечных элементов, N — число всех степеней свободы до задания граничных условий первого рода.
1 Работа поддержана РФФИ (гранты Ж№ 13-01-00101-а, 14-01-00055-а).
На гибридной архитектуре процесс сборки системы уравнений (1.1) возможно осуществить, используя различные вычислительные ресурсы: центральные (CPU), графические (GPU) или одновременно центральные и графические процессоры (CPU+GPU). При выполнении сборки (1.2) на центральном процессоре матрицы Ke, ранее вычисленные на GPU, необходимо переслать в оперативную память. После задания граничных условий первого рода получаем
Ku = f; K = KT; (Kx,x) > 0Ух = 0; K € RNxN,u € RN,f € RN, (1.3)
где N — число неизвестных степеней свободы. Для параллельного решения полученной системы уравнений (1.3) ее необходимо разделить на блоки и разослать по вычислительным ресурсам кластера.
В случае сборки (1.2), выполняемой на нескольких GPU или на CPU+GPU, система уравнений (1.1) сразу формируется в распределенном виде. При этом в соответствующих областях
K
f
необходимо изменить существующие алгоритмы задания условий Дирихле.
В качестве более общего случая будем рассматривать алгоритмы [1] задания граничных условий первого рода на примере задания неоднородного условия Дирихле щ = щ = 0, где i € [1, N — номер неизвестного (степени свободы) в случае распределенной системы (1.1).
В первом алгоритме, использующем умножение на большое число а, полагается кц ^ акц и fi = akaUi, как правило, принимается а ^ 108. В случае распределенной системы (1.1) алгоритм задания граничных условий выполняется на том вычислительном ресурсе, в области памяти которого оказался элемент кц. Следует учесть, что в результате решения (1.3) при использовании первого алгоритма получаем щ ~ щ, так как в строке i матрицы K остаются недиагональные ненулевые элементы. Кроме того, выбор а ограничен сверху обусловленностью системы уравнений (1.3).
Рассмотрим второй алгоритм, в котором над элементами матрицы K и вектоpa f выполняются следующие операции: кц = 1, fi = щ, fj <— f) — kjiUi, kij = 0, kji = 0, j = 1 ,N, здесь N = N. Когда строка i и строки j распределенной матрицы KT, а также соответствующие элементы f оказываются на разных вычислительных узлах, необходимо или переслать ui i
из связей степеней свободы, так, чтобы предотвратить возникновение такого случая. Выбор более эффективного способа реализации зависит от числа степеней свободы, в которых заданы условия Дирихле и архитектуры вычислительной системы.
ui = ui
Уменьшение размерности (N < N) системы уравнений (1.3) требует переупорядочить как уравнения, так и неизвестные, распределенные по областям памяти вычислительной системы. Для снижения затрат на переупорядочение степени свободы нумеруются так, чтобы заданные оказались в конце вектора к или его распределенных блоков. Запишем этот (третий) алгоритм в матричном виде [3]:
Ku = С, K = ctKС, f = ct(f - Ku), (i.4)
здесь K € MNVxNV; u € MNV — вектор искомых степеней свободы; f € MNV — новый вектор правой части; U € MNV — вектор заданных значений степей ей свободы; C € ZxN — матрица перестановок, извлекающая искомые степени свободы; N = N — N&, где Nu — число степеней свободы с заданными условиями Дирихле.
Элементы матрицы С = [Cij ] принимают значения Cj = 1, если Ci = щ-, здесь Ci € u, Uj € щ Cij = 0 — в другом случае. Каждый столбец матрицы С содержит один элемент Cj = 1-Вследствие разреженности матриц С и Ст теоретическая сложность вы числения K в (1.4) составляет O(N-N+N2). Векторы решения систем (1.3) и (1.4) связаны выражением щ = Сщ+щ.
Остановимся на особенностях хранения матриц в приведенных алгоритмах. Второй и третий алгоритмы не предполагают хранения матрицы K без нулевых элементов, потому что в этих случаях требуется удалить часть ненулевых элементов матрицы. Так, во втором алгоритме удаляются ненулевые элементы в строках K, что требует существенных изменений в структуре
данных хранения матрицы, а в третьем алгоритме, кроме того, удаляются строки. Поэтому в этих алгоритмах в компактных форматах хранятся матрицы К, К С и Соответственно, алгоритм с умножением на большое число такого ограничения не имеет.
Применяя первый и второй из рассмотренных алгоритмов к локальным матрицам Ke € € RWeXWe и векторам /е € RWe отдельных конечных элементов е = 1 , m, получим соответствующие поэлементные варианты алгоритмов. Здесь m — число конечных элементов, a Ne — число степеней свободы в конечном элементе е. Задание граничных условий Дирихле на уровне локальных матриц и векторов позволяет применять так называемые matrix-free схемы решения системы (1.3), в которых матрица К не формируется. Например, поэлементные схемы [6], в которых сборка (1.2) заменяется сборкой векторов на этапе решения системы.
§ 2. Учет граничных условий Дирихле в ходе решения системы уравнений
Как правило, существующие варианты задания граничных условий Дирихле в процессе решения конечно-элементной системы линейных алгебраических уравнений (СЛАУ) связаны с матрично-векторными произведениями в явных методах подпространств А.Н. Крылова [4]. В случае неявных методов граничные условия Дирихле должны быть заданы к моменту построения предобуславливателя, а значит до начала решения системы уравнений.
Рассмотрим вариант на основе третьего алгоритма из предыдущего параграфа. Правую часть (1.4) преобразуем до решения системы уравнений. Далее произведение KC вычисляемое на каждой итерации метода решения СЛАУ, запишем в виде последовательности матрично-векторных произведений: р = Ce к = Кр, C = CTq. Таким образом, матричные произведения из (1.4) заменяются на матрично-векторные, которые благодаря разреженности матриц С и
CT
имеют сложность O(N) и O(N) соответственно.
На примере метода сопряженных градиентов (АКТ) рассмотрим более известный вариант задания граничных условий Дирихле [4]. Если задано значение щ = щ = 0 степени свободы г, то полагается, что компоненты г векторов р и q равны нулю на каждой итерации МСГ, здесь q = Кр, р = {pi}, q = {qi}> p,q € В случае неоднородных граничных условий Дирихле столбец г матрицы КС умножается на щ и вычитается из вектора правой части f. Так поступаем с каждым заданным значением степени свободы. Далее, на каждой итерации метода сопряженных градиентов полагаем pi = 0 и qi = 0. Этот подход также применим к поэлементным схемам решения конечно-элементных систем.
Следует отметить, что даже для диагонального предобуславливания требуется сборка диа-
K
алгоритме не будем заменять значения кц на 1, а в правой части зададим fi = щки. В таком варианте diag (К) = diag(K), что позволяет построить и достаточно эффективно применить на гибридной архитектуре диагональный предобуславливатель.
§ 3. Особенности программной реализации
В зависимости от выбранной вычислительной схемы отображение конечно-элементных вычислений по процессорам в том или ином виде наследует разделение расчетной сетки [2]. Для поэлементных схем результатом разделения являются подмножества из m^ г = 1,... ,np, конечных элементов, пересекающихся по общим узлам сетки, здесь np — число параллельных процессов, равное числу центральных процессоров. Это разделение наследуют подмножества матриц конечных элементов {Ke}m=i, вычисление которых распределяется между CPU и GPU. Если численное интегрирование локальных матриц осуществляется на GPU, в том числе на CPU и GPU [5], то матрицы Ke размещаются в памяти ускорителей (все или часть матриц) в массиве для хранения соответствующего подмножества матриц.
Рассмотрим отображение вычислений по параллельным процессам гибридной вычислительной системы. На кластерах гибридной архитектуры реализуется гибридная модель параллелизма (обмен сообщениями + общая память + параллелизм данных), которая, как правило, реализуется при помощи нескольких технологий параллельного программирования: MPI, OpenMP, CUDA. В рамках этой модели каждому подмножеству конечных элементов ставим в соответствие один из np процессов MPI, в котором создается несколько нитей OpenMP (равное числу
ядер CPU). Мастер-нить OpenMP обеспечивает MPI-коммуникации между вычислительными узлами кластера. Другие нити обеспечивают параллелизм при передаче данных между ядрами CPU и несколькими GPU внутри узла, запуск kernel-функций CUDA и вычисления на ядрах центрального процессора.
Кратко остановимся на реализации поэлементного варианта второго алгоритма. На GPU условия Дирихле задаются при помощи kernel-функции CUDA, содержащей программный код вычислений для отдельного конечного элемента. При запуске kernel-функции каждой нити CUDA с номером е ставится в соответствие матрица Ke и вектор fe конечного элемента е, находящегося в памяти ускорителя. Результатом выполнение kernel-функций является множество матриц {Ke}m=1 и вектор /, составленный из векторов fe с заданными условиями Дирихле.
Относительно схем с блочным разделением выше было отмечено, что вычисления матрицы К и вектора / в (1.4) состоят из двух произведений матриц, двух матрично-векторных произведений и разности векторов, все перечисленные операции возможно эффективно выполнить на GPU. Вместе с тем при реализации третьего алгоритма на этапе формирования системы в памяти графического ускорителя необходимо кроме блоков матрицы К хранить также соответствующие блоки матриц К, C и CT, а при решении прикладных задач возможны случаи, когда Nu ^ N и, соответственно, N ~ N. Этот вариант задания условий первого рода более перспективен на этапе решения системы уравнений, где дополнительно потребуется хранить только блоки матриц C и что благодаря разряженности C и CT не существенно отличается от хранения только блоков К и поэтому требует дальнейших исследований.
Как показал опыт исследований, граничные условия первого рода предпочтительнее задавать на этапе формирования конечно-элементных систем. Таким образом, предлагаемые варианты алгоритмов задания граничных условий Дирихле и их реализации на GPU исключают пересылки подмножеств матриц Ke и блоков матрицы К между вычислительными узлами, оперативной памятью и памятью ускорителя, сборку системы уравнений и ее разделение на блоки, увеличивая ускорение конечно-элементных вычислений на гибридных архитектурах в десятки раз, что подтверждается результатами вычислительных экспериментов.
Список литературы
1. Норри Д., де Фриз Ж. Введение в метод конечных элементов. М.: Мир. 1981. 304 с.
2. Копысов С.П., Новиков А.К. Метод декомпозиции для параллельного адаптивного конечно-элементного алгоритма // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2010. № 3. С. 141-154.
3. Даутов Р.З. Программирование МКЭ в MATLAB. Казань: КГУ, 2010. 71 с.
4. Шабров H.H. Метод конечных элементов в расчетах деталей тепловых двигателей. Л.: Машиностроение, 1983. 212 с.
5. Новиков А.К., Копысов С.П., Кузьмин U.M.. Недожогин Н.С. Формирование матриц конечно-элементных систем в GPGPU // Вестник Нижегородского университета им. H.H. Лобачевского. 2014. № 3(1). С. 120-125.
6. Копысов С.П., Кузьмин U.M.. Недожогин Н.С., Новиков А.К., Рычков В.Н., Сагдеева Ю.А., Тон-ков Л.Е. Параллельная реализация конечно-элементных алгоритмов на графических ускорителях в программном комплексе FEStudio // Компьютерные исследования и моделирование. 2014. Т. 6. № 1. С. 79-97.
Поступила в редакцию 01.10.2015
Новиков Александр Константинович, к. ф.-м. н., старший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: [email protected]
Кузьмин Игорь Михайлович, младший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: [email protected]
Недожогин Никита Сергеевич, младший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: [email protected]
A. K. Novikov, I. M. Kuz'min, N. S. Nedozhogin
Several features of the finite element technologies for hybrid computational architecture
Keywords: parallel computation, finite elements method, Dirichlet boundary conditions, hybrid computational architecture.
MSC: 65M60, 65Y05
We consider the features of parallel finite element computation using hybrid architectures related to the boundary conditions of the first kind (the Dirichlet conditions). Variants of the Dirichlet conditions setting in finite element system of equations forming and solving are compared to ensure the preservation of the symmetry of the coefficient matrix.
REFERENCES
1. Norrie D.H., de Vriz G. An introduction to finite element analysis, London: Academic Press, 1978, 301 p. Translated under the title Vvedenie v metod konechnykh elementov, Moscow: Mir, 1981, 304 p.
2. Kopysov S.P., Novikov A.K. Domain decomposition for parallel adaptive finite element algorithm, Vestnik Udmurtskogo Universiteta.. Matematika. Mekhanika. Komp'yuternye nauki, 2010, no. 3, pp. 141-154 (in Russian).
3. Dautov R.Z. Programmirovanie MKE v MATLAB (FEM programming in MATLAB), Kazan: Kazan State University, 2010, 71 p.
4. Shabrov N.N. Metod konechnykh elementov v raschetakh detalei teplovykh dvigatelei (Finite element method in the calculation of parts of heat engines), Leningrad: Mashinostroenie, 1983, 212 p.
5. Novikov A.K., Kopysov S.P., Kuzmin I.M., Nedozhogin N.S. Forming matrices of finite element systems in GPGPU, Vestnik Nizhegorodskogo Universiteta im. N.I. Lobachevskogo, 2014, no. 3(1), pp. 120-125 (in Russian).
6. Kopysov S.P., Kuzmin I.M., Nedozhogin N.S., Novikov A.K., Rychkov V.N., Sagdeeva Yu.A., Tonkov L.E. Parallel implementation of a finite-element algorithms on a graphics accelerator in the software package FEStudio, Komp'yuternye Issledovaniya i Modelirovanie, 2014, vol. 6, no. 1, pp. 79-97 (in Russian).
Received 01.10.2015
Novikov Aleksandr Konstantinovich, Candidate of Physics and Mathematics, Associate Professor, Department of Computational Mechanics, Udmurt State University, ul. Universitetskaya, 1, Izhevsk, 426034, Russia; Senior Researcher, Institute of Mechanics, Ural Branch of RAS, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia.
E-mail: [email protected]
Kuz'min Igor' Mikhailovich, Junior Researcher, Institute of Mechanics, Ural Branch of RAS, ul. T. Baramzi-noi, 34, Izhevsk, 426067, Russia. E-mail: [email protected]
Nedozhogin Nikita Sergeevich, Junior Researcher, Institute of Mechanics, Ural Branch of RAS, ul. T. Baramzi-noi, 34, Izhevsk, 426067, Russia. E-mail: [email protected]
