УДК 004.42, 004.051 Дата подачи статьи: 30.09.22, после доработки: 10.11.22
DOI: 10.15827/0236-235X.141.130-143 2023. Т. 36. № 1. С. 130-143
Векторизация трехмерного метода погруженных границ для повышения эффективности расчетов на микропроцессорах Intel
А.А. Рыбаков 1, к.ф.-м.н., ведущий научный сотрудник, rybakov@jscc.ru А.О. Мещеряков 1, младший научный сотрудник, alex2501@jscc.ru
1 Межведомственный суперкомпьютерный центр РАН (МСЦ РАН) -филиал ФНЦ НИИСИ РАН, г. Москва, 119991, Россия
Работа посвящена повышению эффективности выполнения современных расчетных приложений на высокопроизводительных вычислительных системах. В качестве инструмента повышения эффективности рассматривается векторизация программного кода. С ее помощью однотипные скалярные операции объединяются в векторные аналоги, кратно повышая производительность. Целевой платформой являются современные микропроцессоры Intel, для которых поддержан уникальный набор векторных инструкций AVX-512.
Предлагается подход к векторизации газодинамического решателя, использующего метод погруженных границ и противопотоковую схему Steger-Warming в трехмерном виде. Решатель обладает сложным программным контекстом, автоматическая векторизация которого невозможна. Рассматриваются реализация решателя, а также подходы к организации кода и приведению его к виду, пригодному для автоматической векторизации компилятором icc.
Для обеспечения автоматического применения векторизации к программному коду решателя были применены три основных эквивалентных преобразования. Во-первых, вычисления, одинаковые для всех итераций проведения расчетов, включая матричные операции, были локализованы и вынесены на этап подготовки вычислений. Во-вторых, основные функции решателя были организованы в виде плоских циклов, а структуры данных представлены в виде наборов массивов. В-третьих, к гнездам циклов была применена оптимизация расщепления по условию, с помощью которой можно уменьшить степень разветвленности управления внутри тела цикла. Данные преобразования позволяют компилятору автоматически применять векторизацию кода.
В результате выполненной работы достигнуто ускорение решателя в три раза за счет векторизации при вычислениях на вещественных числах двойной точности.
Ключевые слова: векторизация, оптимизация, газовая динамика, метод погруженных границ, AVX-512.
При численном решении задач газовой динамики часто приходится сталкиваться с расчетными областями со сложной или изменяющейся геометрией (например, обтекание различных тел или протекание потока внутри элементов силовых установок летательных аппаратов). Во время таких расчетов основные проблемы связаны с вычислениями, выполняемыми вблизи сложных границ, - для сложных границ крайне проблематично построение согласованных расчетных сеток. Проблема отчасти может быть решена путем использования неструктурированных и гибридных сеток. В работе [1] представлен подход к генерации гибридной расчетной сетки для численного решения вязкого течения в условиях обтекания тела со сложной геометрией. В основе данного подхода лежит применение триангуляции поверхности обтекаемого тела для построения
расчетной сетки пристеночного слоя, формируемой из слоев согласованных с поверхностью призматических ячеек, тогда как для остального объема используется простая неструктурированная расчетная сетка. Главной проблемой данного подхода является сложность создания слоев призматических ячеек с помощью метода наступающего слоя, при котором необходимо рассчитывать направления движения узлов сетки, размеры шагов и предотвращать самопересечения фронта при обработке вогнутых областей.
Альтернативой является метод погруженных границ [2], с помощью которого можно избежать трудоемкого процесса создания и обслуживания согласованных на поверхности расчетных сеток. Выделим два основных подхода к разрешению граничных условий в методах погруженных границ, различающихся спо-
собом выполнения расчетов на границе: задание граничных условий посредством внешних (источниковых) членов и методы, использующие фиктивные (вспомогательные в расчетах) ячейки (или узлы). При первом из этих подходов в дифференциальные уравнения, описывающие эволюцию течения жидкости или газа, добавляются дополнительные члены, влияние которых подобно присутствию твердого обтекаемого тела, тем самым обозначая его границу. В методе штрафных функций Бринкмана для моделирования обтекаемого твердого тела используется модель пористой среды с приближающимся к нулю значением проницаемости [3].
Вторым подходом к реализации метода погруженных границ является применение так называемых фиктивных расчетных ячеек, которые используются только для вычисления газодинамических величин в других ячейках расчетной сетки. Значения физических величин в фиктивных ячейках вычисляются с помощью аппроксимации по соседним ячейкам расчетной сетки [4]. Как правило, метод погруженных границ с применением фиктивных ячеек подразумевает использование декартовых расчетных сеток, чтобы избежать усложнения шаблонов аппроксимации величин в фиктивных ячейках. Методы, использующие фиктивные ячейки, различаются способами аппроксимации значений в фиктивных ячейках. В работе [5] значения физических величин привязаны к узлам расчетной сетки, а значения в фиктивных узлах, находящихся внутри обтекаемого тела, аппроксимируются с помощью полинома второй степени через так называемую зеркальную точку (то есть точку, симметричную рассматриваемому фиктивному узлу относительно границы раздела среда-тело). В работе [6] рассматривается похожий подход, однако для моделирования течения используются конечно-объемная расчетная сетка и метод Годунова первого порядка с применением приближенного римановского решателя. Таким образом, газодинамические величины привязаны к центрам ячеек, а для аппроксимации значений в центре фиктивной ячейки используются взятые с противоположным знаком значения в зеркальной точке. Другим отличием от метода, предложенного в работе [5], является принудительное отодвижение зеркальной точки от границы в случае, если центр фиктивной ячейки лежит слишком близко к границе, для обеспечения более точной аппроксимации газодинамических величин в ней. Таким обра-
зом, в работах [5, 6] задача аппроксимации значений в фиктивных ячейках сводится к задаче аппроксимации значений внутри поля течения жидкости или газа уже без учета наличия границы.
В настоящей работе рассматривается метод погруженных границ с фиктивными ячейками, используя который, можно ограничиться простыми структурированными расчетными сетками и даже декартовыми сетками с ячейками правильной формы. Основной идеей данного метода является построение структурированной сетки, которая может пересекать сложную границу, при этом некоторые ячейки данной сетки оказываются полностью внутри расчетной области, тогда как другие - только частично. Для ячеек, целиком находящихся в расчетной области, вычисления проводятся привычным образом, а для ячеек, пересекающих границу, выполняются дополнительные вычисления, связанные с учетом граничных условий. Чтобы осуществить учет граничных условий в ячейках, пересекающих расчетную область, необходимо провести коррекцию газодинамических параметров в этих ячейках на основании данных соседних ячеек - выполнить аппроксимацию величин, участвующих в вычислениях. Использование метода погруженной границы с фиктивными ячейками позволяет упростить логику программы, перестроить программный код таким образом, чтобы он оказался более пригодным для декомпозиции задачи, а также для повышения производительности вычислений путем векторизации. За исключением подготовительных действий, а также действий, связанных с аппроксимацией газодинамических параметров вблизи сложной границы, вычисления проводятся таким образом, словно работа ведется просто с обычной структурированной расчетной сеткой. В отличие от рассмотренных методов, использующих зеркальные точки, предлагаемый в данной работе метод выполняет аппроксимацию величин в центрах фиктивных ячеек непосредственно через значения в окружающих ячейках и систему уравнений, порождаемую конфигурацией аппроксимационного шаблона и граничными условиями.
Описание трехмерного метода погруженной границы с использованием фиктивных ячеек
Реализация метода погруженной границы с фиктивными ячейками начинается с построе-
ния структурированной расчетной сетки и разметки ячеек, которые будут участвовать в вычислительном процессе. Для расчетов будем использовать наиболее простой вид сеток - де-картову сетку с кубическими ячейками. Разметка ячеек производится в две фазы. На первой фазе ячейки делятся на четыре типа: обычные (COMMON) - ячейки, полностью расположенные в расчетной области, в них выполняется расчет газодинамических параметров; внутренние (INNER) - ячейки, полностью расположенные внутри обтекаемого тела, в них не происходит никаких газодинамических процессов; граничные (BORDER) - ячейки, пересекающие сложную границу, однако большей частью находящиеся вне обтекаемого тела; фиктивные (GHOST) - ячейки, пересекающие сложную границу, большая часть которых все же находится внутри обтекаемого тела (рис. 1, слева). Условно пригодными для расчетов будем считать ячейки COMMON и BORDER, тогда как ячейки GHOST и INNER относятся, скорее, к обтекаемому телу. Разметка ячеек является самостоятельной задачей и сводится к задаче пересечения объемной и поверхностной расчетных сеток, которая также может быть эффективно векторизована [7], однако в данной работе этот аспект реализации метода не рассматривается. Для разграничения ячеек типов BORDER и GHOST вместо оценки доли объема ячейки, лежащей внутри или вне обтекаемого тела, будем также использовать более простой критерий: если центр ячейки лежит вне тела, она помечается как BORDER, в противном случае - как GHOST (данное упроще-
ние уместно при условии, что берется достаточно мелкая расчетная сетка, линейный размер ячеек в которой меньше изломов поверхности).
После первой фазы разметки ячеек требуется провести корректировку. Во время выполнения расчетов конечно-объемным методом для ячеек COMMON и BORDER нужно будет вычислять потоки вещества, импульса и энергии через все границы этих ячеек, а значит, все их соседи должны содержать газодинамические параметры. Поэтому на фазе коррекции расчетных ячеек все соседи ячеек COMMON и BORDER, которые размечены как INNER, должны быть переквалифицированы в GHOST. После выполнения этого действия нет необходимости различать ячейки COMMON и BORDER, а значит, все ячейки BORDER можно переквалифицировать в COMMON. Таким образом, после выполнения фазы коррекции значение имеют только два типа ячеек: COMMON, для которых проводятся газодинамические расчеты, и GHOST, для которых газодинамические расчеты не проводятся, а выполняется аппроксимация газодинамических параметров с учетом значений соседних ячеек COMMON и граничных условий (рис 1, справа). На рисунке 2 продемонстрирован результат разметки ячеек для тестовой расчетной области (прямоугольная расчетная область, внутри которой выполняется расчет процесса обтекания случайно расположенных сферических объектов различных радиусов).
Центральной задачей реализации метода погруженной границы является выполнение
I Свободный поток Свободный поток
COMMON BORDER BORDER COMMON COMMON COMMON
—
COMMON GHOST INNER COMMON GHOST GHOST
GHOST INNER COMMON GHOST INNER
—' Обтекаемое тело Обтекаемое тело
Рис. 1. Двумерная иллюстрация разметки расчетных ячеек Fig. 1. A two-dimensional illustration of the calculation cell marking
аппроксимации газодинамических величин ОНОБТ-ячеек на основе данных соседних СОММОК-ячеек, а также граничных условий (в данном случае подразумевается граничное условие Неймана). При этом аппроксимацию требуется выполнять как для скалярных величин (плотность и давление), так и для векторной величины скорости. Рассмотрим формулы аппроксимации газодинамических величин для ОИОБТ-ячеек. Данные формулы были получены путем обобщения двумерных уравнений на трехмерный случай, их подробный вывод и описание приведены в [8].
Пусть в пространстве задана скалярная величина ф(х, у, 2). Если известно значение этой величины в четырех точках общего положения фг = ф(хг-, уг, 2г) для I = 0, ..., 3, то можно выполнить ее аппроксимацию в произвольной точке пространства с помощью линейной функции ф(х,у, 2) = ([1, х,у, 2], а), а = [ао, а\, аг, аз]т путем решения системы линейных уравнений:
а0 + а Х0 + а2 Уо + аз 20 =Фо>
а+ал + ал+= ф ,
<
а+ах + а2У2 + а^ = Ф ,
а + аХ + а2Уз + а323 = ф,
то есть а = В^з)ф0123 , где Ф0123 = [ф0, ФЬ Ф2, фз]т, а матрица Б<0ш) имеет следующий вид:
Уо У1 У 2 Уз
При реализации метода погруженных границ для аппроксимации величины ф используем три СОММОК-ячейки (в качестве опорных точек взяты их центры [х1, у1, 21], [х2, у2, 22], [х3, у3, 23]), в которых известны значения этой
величины. В качестве дополнительной точки рассмотрим точку границы [хо, уо, 2о], ближайшую к GHOST-ячейке, а для этой точки границы - граничное условие Неймана, то есть условие равенства нулю производной аппроксимируемой величины по направлению нормали к поверхности в точке [х0, у0, 20]: 5ф , ч , г. 5ф 5ф 5ф
— (Х0: Уо ' 2о) =Ф = 0 = ~ е0, х +— е0, х + ~ е0,х =
de.
dx
dx
dx
= a e + a e + a e
"lc0,x ^ "2C0,y ^ "3C0,z •
Таким образом, вектор коэффициентов аппроксимации определяется как а = В",123 ^0,123 у
0 e0, x e0, У e0,z
1 x1 У1 z1
1 x2 У2 Z2
1 x3 Уз z3
где В 0'12з)
При аппроксимации векторной величины скорости V = [Ух, Уу, Уг] будем считать, что значение этой величины известно в трех соседних СОММОК-ячейках (значения скорости в их центрах равны V , V и V ), а в ближайшей к рассматриваемой GHOST-ячейке точке поверхности (точка с индексом 0) выполняется условие следующего вида: нормальная составляющая вектора скорости и производная тангенциальной составляющей по направлению нормали в точке поверхности равны нулю. Равенство нулю нормальной составляющей V можно записать в виде У0,хе0,х + У0,уе0,у + У0,2^0,2 = 0, откуда, аппроксимировав составляющие вектора V через центр рассматриваемой GHOST-ячейки и центры трех соседних с ней СОММОК-ячеек, получим следующее уравнение относительно неизвестных Уо,х, Уа,у и Уо,г:
(
[1, xo , y0, z0 ]B(G123)
VG, x ^G, У
V1,x \у
e0,x + e0
V 2,x V2, У
V3,x V3,y
после ввода
yG, z
V
V 2, z
V
= 0,
в котором после ввода обозначения
Л = [ЛО, ^2 , ¿з] = [1, Х0 ' У0' 20 ]В(О123) получим
+ ^, уе0,у + ^, уе0,у = "(Е 3=1 ^ Ф ' е0) ) ' ¿О = 2
Что касается условия на тангенциальную составляющую скорости в точке поверхности, то вместо нее будем рассматривать три вектора, перпендикулярные нормали:
7 = [— е0,Уе0,л0] ,
/2 = ["е0,20 е0,Л ,
7 = [0, — е0,2е0,уУ ,
и линеиную аппроксимацию скалярных величин (/, ) с учетом граничного условия Неймана:
е,у + Уе0,х = Т[0, у + \Уе0,х , У2,^0,у +
+У2,уе0,х , -У3,хе0,у + У3,уе0,х ]Т = Тху ,
-УС,Л,г + УС,ге0,х = *е0,г + гб0,х ' -У2,Л,г +
+У2,ге0,х ' -У3,+ У3,ге0,х ]Т = Тхг '
-УС,уе0,г + УС,ге0,у = 1 [0, -У1,уе0,г + \гв0,у ' -У2,уб0,г + +У2,ге0,у ' -У3,уе0,г + У3,ге0,у ]Т = Туг ' где Т = [1, хс > ус > 2С ]В(0'123) •
Таким образом, для определения составляющих скорости сформировалась следующая система уравнений:
"е
yG, y
o" e0 ,z
y o" 0
-e0, z 0
0 -e0, z
Эта система из четырех уравнений содержит три неизвестных, то есть одно уравнение является избыточным. На самом деле избыточным является одно из трех последних уравнений, и из них можно удалить то, в которое входит минимальная по модулю составляющая вектора ё0 .
Несмотря на то, что аппроксимацию газодинамических величин во время расчетов нужно выполнять для каждой GHOST-ячейки на каждой итерации, полученные в приведенных формулах значения векторов d и t (а также матрица из компонентов вектора ё0) зависят только от координат центров ячеек и направлений нормалей к поверхности, то есть изменяются только при изменении геометрии обтекаемого объекта, и в случае неизменяемой геометрии данные значения считаются только один раз на этапе подготовки вычислений. Это существенно, так как только при вычислении этих величин необходимо получать обратные матрицы, что является достаточно дорогой с точки зрения ресурсов процедурой.
Во время выполнения аппроксимации газодинамических величин в GHOST-ячейке возникает вопрос корректности данных вычислений. Аппроксимация может быть выполнена, если существуют обратные матрицы В-^щ и
. Эти матрицы зависят только от выбора
соседних COMMON-точек, по которым выполняется аппроксимация, а также от направления нормали к поверхности обтекаемого тела в
точке поверхности, ближайшей к центру рассматриваемой GHOST-ячейки. Выбранные COMMON-точки для выполнения аппроксимации назовем шаблоном. Существование обратной матрицы В-1щ означает, что точки с индексами G, 1, 2 и 3 должны быть точками общего положения, то есть не лежать в одной плоскости. Существование же обратной матрицы Врщ означает, что нормаль из точки поверхности, ближайшей к рассматриваемой GHOST-ячейке, не должна быть параллельной плоскости, образованной точками с индексами 1, 2 и 3 (которые не лежат на одной прямой).
На рисунке 3 проиллюстрирована схема выбора шаблона аппроксимации газодинамических величин в GHOST-ячейке (центр обозначен черной точкой). Если шаблон выбирается только среди ближайших соседей (индекс ячеек вдоль каждого из направлений структурированной сетки i, j, k отличается не более чем на единицу), то всего существует 26 ячеек, которые претендуют на включение в шаблон. В качестве претендентов рассматриваются только ячейки COMMON, что сразу снижает количество претендентов. Среди оставшихся ячеек нужно перебрать все комбинации троек ячеек и выбрать любую тройку, для которой существуют обратные матрицы В^вщ и .
С большой вероятностью среди ближайших соседей можно найти подходящий шаблон аппроксимации. Если же среди них не удалось найти тройку COMMON-ячеек, которые удовлетворяют всем условиям, то можно расширить поиск.
Рис. 3. Пример корректного (зеленый) и некорректного (красный) шаблонов относительно матрицы
Fig. 3. The example of correct (green) and incorrect (red) patterns relative to the matrix В—
Постановка газодинамической задачи и локализация пригодных для векторизации участков программного кода
Рассмотрим систему уравнений газовой динамики в трехмерном виде:
б б б б — и + — F + — О + — Н = 0, бе бх ду бг
U =
f рм ï
рм рм2 + p
pv , F = pмv
p w p uw
v E У и( E + p) у
G =
E = р
( pv Л ' pw Л
р uv p uw
2 pv + p , H = pvw
pvw 2 pw + p
v v( E + p) y V w(E + p) y
^ V 2
V— + e I, V2 = u2
е (р, р) =-.
( РР ) р(у-1)
Данную систему будем численно решать с помощью метода конечных объемов путем вычисления потоков вектора Г через грани ячейки, перпендикулярные направлению сетки /, вектора О через грани ячейки, перпендикулярные направлению сеткиу, вектора Н через грани ячейки, перпендикулярные направлению сетки к (при вычислении используем декартову сетку с кубическими ячейками, то есть дискретизация пространства выполняется с одним и тем же шагом по всем трем направлениям):
ип+1 = ип- — (1 - F +
и 1 и I (1 I+1/2 1 I-1/2 +
Ah
+Gj+1/2 — Gj-1/2 + Hk +1/2 — Hk-1/2
Для вычисления потоков используем проти-вопоточную схему Steger-Warming [9], в которой потоки получаются следующим образом (на примере потока F):
F+1/2 = F + (un)+f - U",),
f Xj + 2(y - 1)Xj+Xj Л
(u - a)Xj + 2(y - 1)uX j + (u + a)Xj v (Xj + 2(y - 1)Xj+Xj) w(Xj + 2(y - 1)Xj+Xj) (H - ua)Xj + (y - 1)F2Xj + (H + ua)Xj
F=
_P_
2y
= (X( +|X( |) / 2, Xj = и - a, X2 = X3 = X4 = I X5 = и + a, H = (E + p) / р.
С учетом приведенных формул один шаг выполнения расчетов с использованием метода погруженных границ должен содержать следующие необходимые действия.
1. Для всех GHOST-ячеек сетки аппроксимация примитивных газодинамических величин D = [р, u, v, w, p] на основе единожды вычисленных шаблонов (подготовка шаблонов и предварительное вычисление матриц не рассматриваются, так как это однократные действия, вклад которых в общий процесс расчета стремится к нулю при увеличении времени расчетов). Данное действие обозначим approxi-mate_values.
2. Для всех GHOST- и COMMON-ячеек сетки получение из вектора примитивных газодинамических величин D = [р, u, v, w, p] вектора консервативных величин U = [р, pu, pv, pw, E]. Данное действие обозначим d_to_u.
3. Для всех GHOST- и COMMON-ячеек сетки вычисление векторов потоков F±, G± и И±. Эти действия обозначим как calc_fgh.
4. Для всех COMMON-ячеек корректировка вектора консервативных величин U = [р, pu, pv, pw, E] с помощью потоков. Обозначим это действие calc_flows.
5. Для всех COMMON-ячеек обратный пересчет вектора консервативных величин U = [p, pu, pv, pw, E] в вектор примитивных величин D = [p, u, v, w, p]. Данное действие обозначим u_to_d.
Схема выполнения расчетов приведена на рисунке 4.
Использование метода погруженных границ с фиктивными ячейками позволяет проводить газодинамические расчеты вблизи сложных границ без необходимости построения согласованных сеток. С помощью описанного метода и с применением схемы расчета Steger-Warming были проведены расчеты модельной задачи обтекания свободным потоком множества случайных сфер (при этом граница обтекаемых объектов для простоты задавалась аналитически, а не с помощью поверхностной сетки). Реализация модельной задачи находится в публичном репозитории (URL: https://github.com/r-aax/ibm_vec). Для расчетов использовалась простая декартова расчетная сетка с кубическими ячейками. Во время расчетов геометрическая конфигурация оставалась неизменной, что сделало возможным однократное вычисление вспомогательных объектов для аппроксимации газодинамических величин в фиктивных ячейках (в том числе
COMMON- ячейки GHOST- ячейки
Рис. 4. Схема выполнения шага расчетов методом погруженных границ с использованием схемы Steger-Warming
Fig. 4. A scheme for performing a calculation step by the immersed boundary method using the Steger-Warming scheme
речь идет о матрицах В-1щ и ). С помо-
щью метода погруженных границ на декартовой расчетной сетке можно получить качественно правильную картину обтекания, что проиллюстрировано на рисунке 5.
Векторизация и результаты эксперимента
Вычисления на достаточно подробных расчетных сетках с высоким уровнем дискретизации по пространству и времени требуют серьезных вычислительных ресурсов. Для обработки крупных расчетных областей целесообразно применение высокопроизводительных вычислительных кластеров с использованием различных средств повышения быстродействия приложений и распараллеливания вычислений. Основными уровнями распараллеливания являются организация вычислений на распределенной памяти (например, с примене-
нием MPI [10]), использование многопоточного программирования (например, с применением средств OpenMP [11]), а также векторизация программного кода [12], что является низкоуровневой оптимизацией, способной кратно увеличить эффективность исполняемого кода.
Многие современные микропроцессорные архитектуры содержат расширения инструкций для векторизации вычислений. Для ARM наиболее актуальным набором векторных команд на текущий момент является расширенный набор команд под названием NEON -128-битный расширенный SIMD, представленный а архитектурах ARMv7 и ARMv8. Набор команд NEON поддерживает векторные операции по работе с 8-, 16-, 32- и 64-битными знаковыми и беззнаковыми целыми числами, а также с 16-, 32- и 64-битными вещественными числами с плавающей точкой [13]. В современных микропроцессорах архитектуры «Эльбрус» существует поддержка коротких вектор-
Рис. 5. Визуализация поля скоростей (двумерный срез), полученного с помощью метода погруженной границы на расчетной сетке с рисунка 2
Fig. 5. Visualization of the velocity field (two-dimensional slice) obtained using the immersed boundary method on the calculation grid from Fig. 2
ных инструкций, работающих с 8-байтными векторами [14]. Однако наиболее мощная поддержка векторных вычислений реализована в современных линейках микропроцессоров Intel. Это набор векторных инструкций AVX-512, представляющий собой 512-битное расширение 256-битных AVX-инструкций. Набор векторных инструкций AVX-512 содержит огромное количество команд, включая арифметические операции, операции множественного обращения в память, специализированные операции для работы с нейросетями и AES-шиф-рованием, специальные битовые операции, комбинированные операции вида ±a-b±c и многие другие. Однако наиболее значимым преимуществом набора команд AVX-512 является наличие специальных масочных регистров, с помощью которых можно определить подмножество элементов векторов, над которыми требуется произвести векторную операцию. Использование таких масочных регистров позволяет легко векторизовать программный код, представленный в предикатной форме, и, таким образом, выполнять векторизацию кода, содержащего сложное разветвленное управление.
Автоматическая векторизация средствами оптимизирующего компилятора не может быть
применена к программному контексту произвольного вида. Для принудительного применения векторизации могут быть использованы специальные функции-интринсики, которые при компиляции кода раскрываются в конкретные векторные инструкции, последовательности векторных инструкций или библиотечные вызовы (функции-интринсики подключаются к программному коду на языке C/C++ с помощью заголовочного файла <immintrin.h>). Однако прямое ручное создание векторного кода требует серьезных усилий разработчика и такой код сложен для сопровождения. Поэтому потенциальным направлением для повышения эффективности программного кода средствами векторизации является правильная организация кода, чтобы оптимизирующий компилятор мог выполнить векторизацию в автоматическом режиме. Наиболее удобным для выполнения автоматической векторизации программным контекстом является так называемый плоский цикл. Плоским циклом, реализованным на языке программирования C, будем называть цикл for, отвечающий следующим требованиям: все итерации плоского цикла независимы друг от друга (то есть могут выполняться в любом порядке, в том числе и одновременно); внутри плоского цикла осуществляется работа с массивами атомарных данных, причем на i-й итерации цикла к массивам выполняются только обращения вида a\i] (во многих случаях это условие можно выполнить в полном объеме, в других несколько ослабить, жертвуя эффективностью векторизации). Также к программному коду могут предъявляться дополнительные требования, например, выровненность массивов данных в памяти на одну и ту же величину, равную ширине векторного регистра (однако это требование также можно обойти). При организации программного кода должным образом в виде набора плоских циклов компилятор способен выполнить автоматическую векторизацию кода.
В рамках данной работы проводился эксперимент по организации программного кода газодинамического решателя, работающего с данными в формате double (вещественные числа размером 64 бита) с целью векторизации под микропроцессоры Intel Xeon Phi Knights Landing [15]. Данный микропроцессор был выбран, так как, кроме обычных векторных арифметических операций набора AVX-512F, он поддерживает также такие полезные наборы, как AVX-512PF - операции предварительной подкачки данных из памяти и AVX-512ER -
векторные команды для вычисления экспоненты и обратных значений, которые крайне полезны в физических приложениях.
Рассмотрим на примере реализации функции d_to_u из общей схемы вычислений, представленной на рисунке 4, способы организации вычислений интуитивным способом и в виде плоского цикла. Итак, функция d_to_u переводит вектор примитивных газодинамических величин Б = [р, и, V, р] в вектор консервативных величин и = [р, ри, рv, рw, Е] для каждой расчетной ячейки. Интуитивным способом вычисления могут быть организованы в виде следующей схемы: для каждой ячейки создается структура, которая содержит все необходимые величины, и действия по переводу вектора Б в вектор и выполняются над полями данной структуры (листинг 1, слева). При попытке векторизации данного цикла во время объединения нескольких последовательных итераций (объединяться должны 8 итераций, так как размер элемента данных равен 64 битам, а длина векторного регистра 512 битам) вместо чтения из памяти скалярного значения (например,
плотности) должно происходить чтение значений плотности из нескольких последовательно расположенных структур данных. То есть чтение вектора должно осуществляться не из последовательной области памяти. Для таких целей в наборе векторных инструкций AVX-512 предусмотрены инструкции gather/scatter, однако их эффективность гораздо ниже чтения последовательной области памяти размером 512 бит даже при использовании предварительной подкачки данных.
Естественным решением оптимизации вычислений является организация расположения данных в памяти в виде набора массивов (листинг 1, справа). После объединения нескольких итераций цикла команды скалярного доступа в память трансформируются в векторные аналоги доступа к последовательной области памяти размером 512 бит. Вся же арифметика, которая присутствует в коде функции (операции сложения, вычитания, умножения, деления, взятие обратной величины), имеет свои векторные аналоги в наборе инструкций AVX-512.
struct Cell
{ double rho; double rhos[N];
double u; double us[N];
double v; double vs[N];
double w; // организация double ws[N]; // организация
double p; // данных в виде double ps[N]; // данных в виде
double rho u. // массива структур double rho us[N]; // набора массивов
double rho v, double rho vs[N];
double rho w, double rho ws[N];
double E; double Es[N];
Cell cells[N];
void d to u() void d _to_u()
{ for (int i = { 0; i < N; i++) for (int i = 0; i < N; i++)
{ { double rho = cells[i].rho; double rho = rhos[i];
double u = cells[i].u; double u = us[i];
double v = cells[i].v; double v = vs[i];
double w = cells[i].w; double w = ws[i];
double p = cells[i].p; double p = ps[i];
cells[i] rho u = rho * u; rho us[i] = rho * u;
cells[i] rho v = rho * v; rho vs[i] = rho * v;
cells[i] rho w = rho * w; rho ws[i] = rho * w;
cells[i] E = 0.5 * rho Es[i] = 0 .5 * rho
* (u * u + v * v + w * w) * (u * u + v * v + w * w)
} + p / (GAMMA - 1.0); } + p / (GAMMA - 1.0);
Листинг 1. Организация данных для программной реализации функции й (в и в виде массива
структур (слева) и набора отдельных массивов (справа)
Listing 1. Data organization for the software implementation of the d to u function in the form
of an structure array (left) and as a set of separate arrays (right)
Другой критической проблемой, влияющей на эффективность векторизации кода, является наличие условных операций внутри векторизуемого плоского цикла. Конечно, набор инструкций AVX-512 содержит специальные масочные аргументы, с помощью которых можно векторизовать программный код с управлением практически любой сложности (разветвленное управление, гнезда циклов, циклы с вызовами функций, циклы с нерегулярным числом итераций) [16], однако при увеличении количества условий в коде эффективность векторизации снижается.
Как пример такого негативного эффекта можно рассмотреть логику гнезда циклов из функции calc_flows, в которой корректируются консервативные величины с помощью потоков через все грани ячейки. Пусть рассматриваемая расчетная область представлена структурированной сеткой размера ИХ х ИУ х N2 по направлениям г, у, к соответственно. Тогда при учете потоков через левую грань каждой ячейки необходимо отдельно обрабатывать случай г = 0, что соответствует граничному условию расчетной области, причем данное граничное условие может быть разных видов (свободное протекание, условие жесткой стенки) (см. листинг 2, слева). Аналогично нужно рассматривать особые случаи для всех шести граней ячейки, что резко увеличивает количество условий внутри цикла и снижает эффективность векторизации.
Можно заметить, что условия обработки границ расчетной области не уникальны для тройки координат ячейки (г, у, к), но являются константными для некоторого среза ячеек сетки. Это означает, что для гнезда циклов такое условие является частично константным, и гнездо может быть разбито по этому условию. На листинге 2 справа продемонстрировано разбиение гнезда циклов по условию г = 0, после чего итоговое гнездо, содержащее основную часть вычислений, освобождается от условия. Аналогичным образом можно выполнить разбиение по остальным условиям, освободив от них основное гнездо циклов, которое после этого успешно векторизуется.
На рисунке 4 приведена общая схема выполнения расчетов для одной итерации численного метода. Основными функциями в этих расчетах являются approximate_values, d_to_u, са1с_£§^ calc_flows, u_to_d. Во время проведения эксперимента на модельной задаче был собран профиль исполнения до векторизации и после нее, результаты распределения времени исполнения между этими основными функциями представлены на рисунке 6.
Результаты ускорения после векторизации каждой функции в отдельности и суммарного ускорения всего расчетного кода отражены на рисунке 7. Заметим, что во время векторизации расчетного кода, работающего с вещественными числами двойной точности, один векторный регистр содержит 8 элементов данных.
approximate_values Bd_to_u calcjgh ■ calcjlows u_to_d ■ approximate_values ■ d_to_u calcjgh Bcalcjlows ■ u_to_d
Рис. 6. Распределение времени выполнения между основными функциями до (слева)
и после (справа) векторизации
Fig. 6. Distribution of execution time between the main functions before (left) and after (right) vectorization
4 3 2 1 0
approximate_values
d to u
calc_fgh
calc flows
u_to_d
Рис. 7. Ускорение отдельных функций и суммарное ускорение Fig. 7. Acceleration of individual functions and total acceleration
6
5
Можно эмпирически считать, что при замене в плоском цикле всех скалярных инструкций на точные векторные аналоги логическое ускорение составит 8 раз (количество инструкций в итерации цикла останется неизменным, но количество итераций в плоском цикле уменьшится в 8 раз). Однако на практике ускорение в 4 раза уже можно считать успешным.
Следует заметить, что ускорения различных функций в результате векторизации сильно отличаются друг от друга. Причиной являются конкретные особенности каждой из этих функций.
Наибольшее ускорение продемонстрировала функция calc_fgh, реализующая вычисление потоков Е±, С±, И±. Этому способствовали следующие ее особенности. Во-первых, функция содержит плоский цикл, полностью удовлетворяющий всем условиям, что делает возможной замену всех скалярных операций на векторные аналоги. Во-вторых, в скалярной
версии используются библиотечные вызовы abs (вычисление абсолютного значения) и sqrt (вычисление квадратного корня). В наборе векторных инструкций AVX-512 есть векторные операции VAND (с помощью побитовой операции над вещественным числом можно обратить в ноль бит знака данного числа) и VSQRT, которые реализуют эти действия просто одной операцией, что сильно ускоряет исполнение. Наконец, обилие операций умножения и сложения делает возможным применение векторных комбинированных операций вида ±а -Ь±е.
Ускорение выше среднего продемонстрировали функции d_to_u и u_to_d. Они также состоят из плоских циклов, поэтому скалярные операции могут быть заменены на векторные аналоги. Более низкое ускорение функции u_to_d объясняется наличием операций деления, которые выполняются медленнее операций сложения и умножения.
Показатель ускорения функции calc_flows составил менее двух раз даже после избавления от всех условий, связанных с обработкой границ расчетной области. Такая низкая эффективность векторизации объясняется тем, что внутри этой функции присутствуют циклы, которые не являются плоскими. Например, в функции calc_flows для корректировки консервативных величин внутри ячейки (i, j, к) требуется обращаться за данными ко всем смежным по граням ячейкам: (i ± 1, j, к), (i, j ± 1, k), (i, j, к ± 1). Это нарушает требование по унификации обращения к массивам данных на одной итерации плоского цикла и немедленно приводит к появлению операций gather/scatter, понижая эффективность векторизации. Для функции approximate_values ситуация выглядит аналогичной, так как для выполнения аппроксимации данных в ячейке (i, j, k) требуется обращаться за данными в ячейки, относящиеся к шаблону аппроксимации (на рисунке 4 нарушение требования унификации обращения за данными продемонстрировано красными стрелками, означающими, что при обработке одной ячейки сетки возникает необходимость обращаться за данными к другой ячейке).
Несмотря на низкие показатели эффективности векторизации функций, содержащих циклы, не являющиеся плоскими, общая эффективность от векторизации метода погруженных границ составила примерно три раза, при этом удалось ограничиться только реорга-
низацией кода без использования ассемблера или функций-интринсиков для форсированного использования векторных инструкций.
Заключение
В работе рассмотрен метод погруженных границ для выполнения газодинамических расчетов обтекания тел со сложной геометрией. Данный метод реализован в трехмерной постановке в условиях постоянной геометрии и аналитически заданных границ обтекаемых тел, разработан программный код на языке C. Для тестовой задачи проведены расчеты с использованием противопотоковой схемы Steger-Warming.
Изучены вопросы векторизации программного кода для микропроцессоров Intel с поддержкой набора инструкций AVX-512, предложены подходы по реорганизации программного кода для обеспечения автоматической векторизации кода с помощью оптимизирующего компилятора icc для микропроцессора Intel Xeon Phi Knights Landing. В результате экспериментов достигнуто ускорение расчетного кода за счет векторизации при использовании вещественных чисел двойной точности.
Предложенные методы подготовки кода для векторизации могут быть использованы для повышения эффективности любых расчетных приложений при представлении логики программы в виде объединения плоских циклов.
Работа выполнена в МСЦ РАН в рамках госзадания по теме FNEF-2022-0016. В исследованиях использовался суперкомпьютер МВС-10П.
Литература
1. Ye H., Liu Y., Chen B., Liu Z., Zheng J., Pang Y., Chen J. Hybrid grid generation for viscous flow simulations in complex geometries. Advances in Aerodynamics, 2020, vol. 2, art. 17. DOI: 10.1186/s42774 020-00042-x.
2. Цветкова В.О., Абалакин И.В., Бобков В.Г., Жданова Н.С., Козубская Т.К., Кудрявцева Л.Н. Моделирование обтекания винта на адаптивной неструктурированной сетке с использованием метода погруженных границ // Математическое моделирование. 2021. Т. 33. № 8. С. 59-82. DOI: 10.20948/mm-2021-08-04.
3. Жданова Н.С., Абалакин И.В., Васильев О.В. Расширение метода штрафных функций Бринк-мана для сжимаемых течений вокруг подвижных твердых тел // Математическое моделирование. 2022. Т. 34. № 2. С. 41-57. DOI: 10.20948/mm-2022-02-04.
4. Chi C., Abdelsamie A., Thevenin D. A directional ghost-cell immersed boundary method for incompressible flows. J. of Comput. Phys., 2020, vol. 404, art. 109122. DOI: 10.1016/j.jcp.2019.109122.
5. Yousefzadeh M., Battiato I. High-order ghost-cell immersed boundary method for generalized boundary conditions. Int. J. of Heat and Mass Transfer, 2019, vol. 137, pp. 585-598. DOI: 10.1016/j.ijheatmasstrans-fer.2019.03.061.
6. Chi C., Lee B.J., Im H.G. An improved ghost-cell immersed boundary method for compressible flow simulations. Int. J. for Numerical Methods in Fluids, 2017, vol. 83, pp. 132-148. DOI: 10.1002/fld.4262.
7. Рыбаков А.А. Векторизация нахождения пересечения объемной и поверхностной расчетных сеток для микропроцессоров с поддержкой AVX-512 // Тр. НИИСИ РАН. 2019. Т. 9. № 5. С. 5-14.
8. Рыбаков А.А. Метод погруженной границы с использованием фиктивных ячеек в трехмерной постановке // Современные информационные технологии и ИТ-образование. 2020. Т. 16. № 2. С. 321-330.
9. Смирнова Н.С. Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа // Тр. МФТИ. 2018. Т. 10. № 1. С. 122-141.
10. Chunduri S., Parker S., Balaji P., Harms K., Kumaran K. Characterization of MPI usage on a production supercomputer. Proc. Int. Conf. for High Performance Computing, Networking, Storage and Analysis, 2018, pp. 386-400. DOI: 10.1109/SC.2018.00033.
11. Maj eed R., Farrukh R., Riaz O., Ali S., Samad A., Khan M. Parallel implementation of FEM solver for shared memory. Math. Problems in Engineering, 2022, vol. 2022, art. 2375102. DOI: 10.1155/2022/2375102.
12. Tian X., Saito H., Preis S.V., Garcia E.N., Kozhukhov S.S. et al. Effective SIMD vectorization for Intel Xeon Phi coprocessors. Sci. Programming, 2015, art. 269764. DOI: 10.1155/2015/269764.
13. Arm Developer. Introducing NEON Development Article. URL: https://developer.arm.com/documen-tation/dht0002/a/Introducing-NEON/What-is-NEON- (дата обращения: 28.09.2022).
14. Волконский В.Ю., Брегер А.В., Бучнев А.Ю., Грабежной А.В., Ермолицкий А.В., Муханов Л.Е. и др. Методы распараллеливания программ в оптимизирующем компиляторе // Вопросы радиоэлектроники. 2012. Т. 4. № 3. С. 63-88.
15. Jeffers J., Reinders J., Sodani A. Intel Xeon Phi Processor High Performance Programming. Morgan Kaufmann Publ., 2016, 662 p.
16. Savin G.I., Shabanov B.M., Rybakov A.A., Shumilin S.S. Vectorization of flat loops of arbitrary structure using instructions AVX-512. Lobachevskii J. of Math., 2020, vol. 41, no. 12, pp. 2566-2574. DOI: 10.1134/S1995080220120331.
Software & Systems Received 30.09.22, Revised 10.11.22
DOI: 10.15827/0236-235X.141.130-143 2023, vol. 36, no. 1, pp. 130-143
Vectorization of the three-dimensional immersed boundary method for improving the efficiency of calculations on Intel microprocessors
A.A. Rybakov 1, Ph.D. (Physics and Mathematics), Leading Researcher, rybakov@jscc.ru A.O. Meshcheryakov 1, Junior Researcher, alexmeshr@jscc.ru
1 Joint Supercomputer Center of the Russian Academy of Sciences - JSCC, Moscow, 119991, Russian Federation
Abstract. The work is devoted to increasing the efficiency of modern computational applications on highperformance computing systems. The authors consider program code vectorization as a tool for increasing efficiency. Vectorization helps combining scalar operations of the same type into vector analogs, significantly increasing performance. Modern Intel microprocessors were chosen as the target platform, for which a unique set of vector instructions AVX-512 is supported.
The paper considers an approach to vectorization of a gas dynamic solver using the immersed boundary method and the Steger-Warming upwind scheme in 3D case. This solver has a complex programming context that cannot be vectorized automatically. The paper considers the implementation of the solver, as well as approaches to organizing the code and bringing it to a form suitable for automatic vectorization by the icc compiler.
To ensure automatic application of vectorization to the solver code, three basic equivalent transformations were applied. First, the calculations, which are the same for all iterations including matrix operations, were localized and brought to the stage of preparing calculations. Second, the main functions of the solver were organized as flat loops, and the data structures were presented as sets of arrays. Third, splitting by condition optimization was applied to loop nests, which can be used to reduce the degree of control branching inside the loop body. These transformations allow the compiler to automatically apply code vectorization.
As a result of the work performed, the solver was accelerated by a factor of 3 due to vectorization when performing calculations on double-precision real numbers.
Keywords: vectorization, optimization, fluid dynamics, immersed boundary method, AVX-512.
Acknowledgements. The work was carried out at the Joint Supercomputer Center of the Russian Academy of Sciences within the framework of the state assignment on the topic FNEF 2022 0016. The authors used MVS-10P supercomputer in the studies.
References
1. Ye H., Liu Y., Chen B., Liu Z., Zheng J., Pang Y., Chen J. Hybrid grid generation for viscous flow simulations in complex geometries. Advances in Aerodynamics, 2020, vol. 2, art. 17. DOI: 10.1186/s42774 020-00042-x.
2. Cvetkova V.O., Abalakin I.V., Bobkov V.G., Zhdanova N.S., Kozubskaya T.K., Kudryavceva L.N. Simulation of flow near rotating propeller on adaptive unstructured meshes using immersed boundary method. Math. Models Comput. Simul., 2021, vol. 33, no. 8, pp. 59-82. DOI: 10.20948/mm-2021-08-04 (in Russ.).
3. Zhdanova N.S., Abalakin I.V., Vasilev O.V. Generalized Brinkman volume penalization method for compressible flows around moving obstacles. Math. Models Comput. Simul., 2022, vol. 34, no. 2, pp. 41-57. DOI: 10.20948/mm-2022-02-04 (in Russ.).
4. Chi C., Abdelsamie A., Thevenin D. A directional ghost-cell immersed boundary method for incompressible flows. J. of Comput. Phys, 2020, vol. 404, art. 109122. DOI: 10.1016/j.jcp.2019.109122.
5. Yousefzadeh M., Battiato I. High-order ghost-cell immersed boundary method for generalized boundary conditions. Int. J. of Heat and Mass Transfer, 2019, vol. 137, pp. 585-598. DOI: 10.1016/j.ijheatmasstrans-fer.2019.03.061.
6. Chi C., Lee B.J., Im H.G. An improved ghost-cell immersed boundary method for compressible flow simulations. Int. J. for Numerical Methods in Fluids, 2017, vol. 83, pp. 132-148. DOI: 10.1002/fld.4262.
7. Rybakov A.A. Vectorization of finding the intersection of volume grid and surface grid for microprocessors with AVX-512 support. Proc. SRISA RAS, 2019, vol. 9, no. 5, pp. 5-14 (in Russ.).
8. Rybakov A.A. Immersed boundary method using ghost cells in a three-dimensional case. Modern Information Technologies andIT-Education, 2020, vol. 16, no. 2, pp. 321-330 (in Russ.).
9. Smirnova N.S. Comparison of flux splitting schemes for numerical solution of the compressible Euler equations. Proc. MIPT, 2018, vol. 10, no. 1, pp. 122-141 (in Russ.).
10. Chunduri S., Parker S., Balaji P., Harms K., Kumaran K. Characterization of MPI usage on a production supercomputer. Proc. Int. Conf. for High Performance Computing, Networking, Storage and Analysis, 2018, pp. 386-400. DOI: 10.1109/SC.2018.00033.
11. Majeed R., Farrukh R., Riaz O., Ali S., Samad A., Khan M. Parallel implementation of FEM solver for shared memory. Math. Problems in Engineering, 2022, vol. 2022, art. 2375102. DOI: 10.1155/2022/2375102.
12. Tian X., Saito H., Preis S.V., Garcia E.N., Kozhukhov S.S. et al. Effective SIMD vectorization for Intel Xeon Phi coprocessors. Sci. Programming, 2015, art. 269764. DOI: 10.1155/2015/269764.
13.Arm Developer. Introducing NEON Development Article. Available at: https://developer.arm.com/ documentation/dht0002/a/Introducing-NEON/What-is-NEON- (accessed September 28, 2022).
14. Volkonsky V.Yu., Breger A.V., Buchnev A.Yu., Grabezhnoy A.V., Ermolicky A.V., Mukhanov L.E. et al. Program parallelization methods implemented in optimizing compiler. Issues of Radio Electronics, 2012, vol. 4, no. 3, pp. 63-88 (in Russ.).
15. Jeffers J., Reinders J., Sodani A. Intel Xeon Phi Processor High Performance Programming. Morgan Kaufmann Publ., 2016, 662 p.
16. Savin G.I., Shabanov B.M., Rybakov A.A., Shumilin S.S. Vectorization of flat loops of arbitrary structure using instructions AVX-512. Lobachevskii J. of Math., 2020, vol. 41, no. 12, pp. 2566-2574. DOI: 10.1134/S1995080220120331.
Для цитирования
Рыбаков А.А., Мещеряков А.О. Векторизация трехмерного метода погруженных границ для повышения эффективности расчетов на микропроцессорах Intel // Программные продукты и системы. 2023. Т. 36. № 1. С. 130-143. DOI: 10.15827/0236-235X.141.130-143.
For citation
Rybakov A.A., Meshcheryakov A.O. Vectorization of the three-dimensional immersed boundary method for improving the efficiency of calculations on Intel microprocessors. Software & Systems, 2023, vol. 36, no. 1, pp. 130-143 (in Russ.). DOI: 10.15827/0236-235X.141.130-143.