Описывается библиотека операторов для работы с сеточными функциями, определёнными на трёхмерных локально-адаптивных сетках. Библиотека спроектирована так, чтобы детали ее реализации были скрыты от пользователя, что позволяет эффективно реализовать ее на машинах с различной архитектурой (параллельных, гибридных и т.п.). При создании библиотеки ставилось несколько целей: приближение внешнего вида программ к формулам в теоретических работах, относительная простота использования, в том числе при работе на графических платах с архитектурой CUDA, вычислительная эффективность.
вые параллельные архитектуры, например, на графические ускорители (например, CUDA, см. [2]) и новые процессоры Intel Xeon Phi (см. [3]). Существует ряд систем, облегчающий подобный перенос. Среди них можно отметить разработанную в ИПМ им. М.В. Келдыша систему DVM (см. [4]), системы OpenMP (см. [5]) и OpenACC (см. [6]). С использованием этих и им подобных систем можно относительно легко преобразовать существующую последовательную программу, написанную на языках C/C++ или Fortran для работы на графических ускорителях или процессорах Intel Xeon Phi. Для этого в программе перед участками кода (как правило, циклами), которые должны исполняться параллельно, нужно вставить специальные инструкции компилятору, которые позволят ему сгенерировать параллельный код. В языках C и C++ это, как правило, делается с помощью т.н. прагм (#pragma dvm, #pragma omp или #pragma acc). В программах на языке Fortran используются псевдокомментарии (например, !$acc). Если скомпилировать такую программу обычным компилятором, который «не понимает» этих инструкций, то он их просто проигнорирует и сгенерирует последовательный код. Всё это позволяет относительно легко перенести существующую программу на новые архитектуры. Если программа компилируется для гетерогенной архитектуры с отдельным ускорителем (например, с графическим ускорителем, таким, как CUDA или с процессором Intel Xeon Phi), то нужно иметь в виду, что у отдельного ускорителя своя оперативная память, и он может обрабатывать только данные, размещённые в ней. Поэтому компилятор, генерирующий параллельный код, как правило, хранит две копии данных – в памяти основной (host) системы и в памяти ускорителя и при необходимости синхронизует их (копирует в ту или другую сторону). Синхронизацией данных (когда и что копировать) также можно управлять с помощью специальных инструкций (прагм или псевдокомментариев).
Однако, всё это никак не сокращает размеров формул в программах.
Описываемая библиотека стоит на стыке этих двух направлений. С одной стороны, она позволяет использовать сокращённую запись выражений с использованием операторов, аналогичных математическим, а с другой стороны, позволяет путём простой перекомпиляции генерировать код, исполняемый на различных архитектурах. 2. Библиотека gridmath
Для того, чтобы внешний вид уравнений в программах на языке C++ можно было
приблизить к внешнему виду тех же уравнений в математической литературе, написана операторная
библиотека gridmath (см. [
с использованием библиотеки gridmath может быть записано, например, так: u=laplas(mu*v);
В правой части оператора присваивания может стоять выражение любой сложности, результатом этого выражения будет сеточный вычислитель, который запомнит всю цепочку вычислений, не производя их. В момент присваивания вычислителя сеточной функции произойдёт запуск вычислений, и для всех точек сеточной функции в левой части оператора присваивания будут произведены вычисления в соответствие с запомненной цепочкой.
Библиотека gridmath решает ещё одну важную задачу. Как известно, на многих
современных суперкомпьютерах (кластерах) стоят графические ускорители CUDA или ускорители Intel
Xeon Phi. Например, в списке top500 самых мощных суперкомпьютеров (см [8]) за ноябрь 2014
года на первом месте стоит суперкомпьютер Tianhe-2 (Китай - NUDT), основанный на
процессорах Intel Xeon Phi 31S1P, а на втором месте – Titan (США - Cray Inc), основанный на
графических ускорителях NVIDIA K20x. Многие Российские суперкомпьютеры, например, самый
мощный Российский суперкомпьютер «Ломоносов» (см. [
Как показано выше, основной запуск вычислений в библиотеке производится в операторе присваивания вычислителя сеточной функции. То, как этот оператор реализован, полностью скрыто от пользователя и может производиться как последовательно, так и параллельно для точек сеточной функции в левой части. Порядок обхода точек – внутреннее дело библиотеки и пользователь не должен закладываться на тот или иной порядок обхода точек. Всё это даёт библиотеке существенную гибкость. В частности, если компиляция осуществляется компилятором nvcc (CUDA), то, во-первых, все данные сеточных функций располагаются в памяти графического ускорителя, а не host-системы, и, во-вторых, оператор присваивания реализован путём вызова ядра (kernel) в графическом процессоре и осуществляется параллельно. В последовательном случае обход точек в случае трёхмерной регулярной сетки осуществляется в трёх вложенных циклах по трём осям. При работе на CUDA следует обратить внимание на то, что, в отличие от упомянутых во введении систем, нацеленных на распараллеливание циклов (DVM, OpenACC), библиотека на хранит копию данных в памяти основного (host) процессора. Все данные сеточных функций хранятся только в памяти CUDA. Для хранения данных сеточных функций при компиляции «обычным» компилятором используется класс std::vector из стандартной библиотеки C++, а при компиляции с помощью nvcc используется класс thrust::device_vector из CUDA SDK.
Таким образом, при использовании библиотеки gridmath перенос программы на графические ускорители CUDA является относительно простой задачей и не требует от пользователя библиотеки (а им чаще всего является математик, специалист по численным методам) глубокого знания методов программирования для CUDA. В простейшем случае программу достаточно просто перекомпилировать. Немного сложнее получается при использовании MPI. В этом случае надо учитывать, что данные сеточных функций расположены в памяти CUDA, а MPI работает с памятью host-системы, поэтому нужно ещё добавить копирование данных из CUDA в host-систему и обратно.
Перенос последовательной программы, написанной с использованием библиотеки
gridmath, на ускоритель Intel Xeon Phi также является простой задачей. Для этого используется
OpenMP (см. [
Первоначальная версия библиотеки написана для регулярных двух и трёхмерных
прямоугольных сеток. С её помощью успешно перенесён на CUDA многосеточный метод (см. [
Затем встал вопрос о переносе библиотеки на нерегулярные сетки. Основная идея такого переноса состоит в том, чтобы некоторым образом упорядочить ячейки сетки, а затем, при присваивании выражения сеточной функции, делать их одномерный обход, последовательный (в одном цикле) или параллельный (для CUDA или OpenMP). При переносе библиотеки на другие типы сеток по сути, создаётся новая библиотека, а не новая версия существующей библиотеки. В результате создаётся целое семейство библиотек, отличающихся типом сетки, но с общей внутренней идеологией и с многочисленными общими частями.
Важное отличие в случае нерегулярных сеток состоит в том, что сетка явно отделяется от сеточной функции. В простейшем случае в этом не было необходимости, т.к. регулярная прямоугольная сетка устроена очень просто и её вынос в отдельный объект не имеет смысла. В случае же нерегулярной сетки это не так. Есть целый ряд нетривиальных операций с сеткой, таких как определение соседей ячейки, вычисление площадей граней, объёмов ячеек, векторов нормалей к граням, не имеющих непосредственного отношения к данным, которые хранит сеточная функция в каждой ячейке. Поэтому вынос сетки в отдельный самостоятельный объект в этом случае оправдан. Сетка не хранит данных, не имеющих к ней непосредственного отношения, этим занимается сеточная функция. Сеточная функция определяется на сетке, т.е. объектсетка передаётся сеточной функции как параметр создания и может впоследствии использоваться. Сеточная функция по желанию программиста может быть определена в ячейках сетки или в её узлах (вершинах). Сеточная функция может хранить для каждой ячейки сетки одно или несколько скалярных или векторных значений. Например, если нужно для каждой ячейки сетки хранить некий трёхмерный вектор (например, вектор скорости), то нет необходимости заводить три сеточных функции для хранения каждой из координат. Можно завести одну сеточную функцию и хранить в ней для каждой ячейки сетки сразу все три координаты. На одной и той же сетке можно определить любое количество сеточных функций. Одна из них может, например, хранить скорость (все три координаты), другая – давление, третья – плотность, четвёртая – температуру и т.д. Библиотека позволяет хранить и размерные величины, при этом размерность переменной – это параметр времени компиляции и никак не влияет ни на объём данных, занимаемой переменной в оперативной памяти, ни на скорость обработки данных. С другой стороны, использование размерных величин, во-первых, повышает наглядность программы, т.к. уже из определения переменной (её размерности) ясно, что в ней содержится, и, во-вторых, позволяет часть программных ошибок (таких, как сложение и вычитание величин, имеющих разную размерность) отловить уже на стадии компиляции.
В качестве первого шага при переносе библиотеки на произвольные нерегулярные сетки было принято решение перенести библиотеку на локально-адаптивные сетки.
Локально адаптивные сетки формируются на основе регулярных равномерных сеток.
Физический шаг исходной регулярной сетки вдоль каждой из осей постоянный, но вдоль разных
осей может быть разным. Сетка строится самой библиотекой или может быть загружена из
файла. После построения сетку можно сохранить в файл. Формат файла бинарный и
нестандартный. Возможен экспорт сетки в файл формата TecPlot (одно зонный и многозонный), а
также в формат ig и METIS (см [
Сетка строится следующим образом. Внутрь исходной регулярной сетки помещается тело произвольной формы, например, шар или цилиндр так, чтобы был простой алгоритм, позволяющий для любой точки внутри исходной физической области ответить на вопрос, находится ли эта точка внутри тела или снаружи. Затем строятся несколько (например, 5) уровней локально-адаптивной сетки, причём каждый следующий уровень строится на основе предыдущего. Исходная регулярная сетка рассматривается как нулевой уровень. До построения сетки формируются глобальные (общие для всех уровней) массивы физических координат по каждому из трёх направлений. Первоначальное заполнение этих массивов простое – значение по i-му индексу равно шагу по данному направлению, помноженному на этот индекс. Эти массивы динамические – по мере построения сетки они расширяются. При делении ребра пополам может возникнуть новая точка с координатой, которой ещё нет в массиве координат. В этом случае координата с заданным значением добавляется в конец массива координат. Это, в свою очередь, означает, что массивы координат не упорядочены по возрастанию (кроме начала этих массивов), значения координат могут идти в произвольном порядке.
Алгоритм построений следующего (начиная с первого) уровня следующий. Изначально для строящегося уровня список рёбер, граней и ячеек пустой. Рассматриваются все рёбра сетки предыдущего уровня. У каждого ребра есть две вершины. Если оказывается так, что одна из вершин ребра лежит внутри тела, а вторая – снаружи, то это ребро делится пополам, и получившиеся два коротких ребра добавляются к списку рёбер строящегося уровня, а само ребро с предыдущего уровня помечается как поделенное. В описании каждого из объектов сетки (вершины, грани или ячейки) хранится индекс соответствующего объекта со следующего уровня. Если этот индекс равен -1 (недопустимое значение для индекса), то это означает, что данный объект сетки не поделён, иначе в нём хранится индекс первого из объектов, из которых состоит данный объект. После того, как все рёбра рассмотрены (и некоторые из них поделены), рассматриваются все грани сетки предыдущего уровня. Если оказывается, что для очередной грани хотя бы одно из четырёх рёбер было поделено на предыдущем шаге, то оставшиеся не поделёнными рёбра этой грани также делятся пополам, а сама грань делится на 4 маленьких равных грани, которые добавляются к списку граней строящегося уровня, а грань предыдущего уровня помечается как поделенная. Затем аналогичная процедура повторяется для ячеек. Если хотя бы одна из шести граней ячейки была поделена на предыдущем шаге, то оставшиеся не поделёнными грани также делятся (вместе со своими рёбрами), а сама ячейка делится на 8 маленьких равных ячеек, которые добавляются к списку ячеек строящегося уровня, а ячейка предыдущего уровня помечается как поделенная.
После построения очередного уровня в том случае, если номер построенного уровня больше 1, осуществляется сглаживание предыдущих уровней в обратном порядке (вначале сглаживается пред-предыдущий уровень, затем пред-пред-предыдущий и т.д.). Сглаживание делается для того, чтобы с двух сторон одной грани не было ячеек «через уровень», т.е. это должны быть или ячейки одного уровня, или соседних уровней. Алгоритм сглаживания следующий. Просматриваются все ячейки сглаживаемого уровня. Если какая-то ячейка не была поделена, а хотя бы одна из её граней поделена дважды (поделена грань и хотя бы одна из её дочерних граней), то ячейка делится (вместе со всеми своими гранями и рёбрами).
Формат внутреннего представления информации обо всех элементах сетки (рёбрах, гранях, ячейках) одинаковый и содержит следующие поля: int child; int boundary_type; int map[6];
Поле child задаёт индекс первого дочернего объектов следующего уровня или -1.
Поле boundary_type содержит признак граничного элемента. 0 – внутренний элемент, 1 – граничный с граничным условием типа Дирихле, 2 - граничный с граничным условием типа Неймана.
В поле map хранится информация об объектах того же уровня, составляющих данный элемент (вершинах для ребра, рёбер для граней, граней для ячеек). Для рёбер хранятся координаты (глобальные индексы) двух вершин. Для граней гранятся индексы 4 рёбер, окружающих данную грань. Последние два элемента массива map в этом случае не используются. Для ячеек хранятся индексы 6 граней, окружающих данную ячейку.
Для сеточных функций, определённых на локально-адаптивных сетках, определены такие же операции, как и для сеточных функций на регулярных сетках (сложение, вычитание, умножение и пр. со скалярами и вычисляемыми объектами) и также можно создавать операторы. 4. Решение параболических уравнений Рассмотрим параболическое уравнение Будем искать решение этого уравнения в трёхмерной области Ω = [x0 ; x1] × [y0 ; y1] × [z0 ; z1] (прямоугольный параллелепипед) при условии, что на всех гранях расчётной области задан нулевой поток , где – внешняя нормаль к границе расчётной области. Функции κ(r), f(r), a(r) ≥ 0 являются заданными, r=(t,x,y,z) ∈ G = [t0 ; T] × Ω, [t0 ; T] – заданный интервал времени.
Предполагается, что расчётная область может содержать зоны со сложной геометрией, отличной от прямоугольных параллелепипедов. Эти зоны могут иметь различные физические свойства, на границах зон тензор диффузии κ может быть разрывным и его скачок на поверхности разрыва может быть большим. Тогда в окрестности поверхности разрыва решение может иметь большие градиенты, для сеточного воспроизведения которых обычно требуются подробные сетки. Но построение подробной сетки на всей расчётной области приведёт к значительному росту вычислительной сложности задачи. Выборочное измельчение элементов сетки может оказаться компромиссным решением. Локально-адаптивные сетки как раз и помогают решить такую задачу.
Рассмотрим уравнение теплопроводности в следующей форме: Для получения дискретных уравнений используется метод конечных объёмов. Сеточная функция задана в центрах ячеек, и ячейка консервативности совпадает с ячейкой сетки: , Индекс i нумерует сеточную функцию в центре ячейки сетки, – искомое значение функции на верхнем временном слое – площадь боковой грани ячейки, – удельный диффузионный поток через боковую грань ячейки в направлении внешней нормали, Vi – объём i-ой ячейки.
Для параболических уравнений универсальной является явная схема, но она используется крайне редко из-за известного ограничения на шаг по времени. В данном примере используется явно-итерационная схема ЛИ-М, которая специально разработана для решения параболических уравнений (см. [14]).
Данная задача решена ранее Феодоритовой О.Б. без использования библиотеки gridmath. Интересно сравнить по быстродействию старую последовательную версию с последовательной же версией, но с использованием библиотеки, а также с CUDA версией (также с использованием библиотеки). Как уже говорилось выше, CUDA версия программы получается практически простой перекомпиляцией компилятором nvcc. Решалась также распределённая задача с обменом по MPI. Декомпозиция сетки на несколько узлов была сделана Головченко Е.Н. с помощью разработанного ею метода (см. [15]).
Для решения задачи разработан разностный оператор диффузии. На вход этому оператору подаются значения тензора диффузии axx, ayy, azz и функции a во всех ячейках сетки. На основе этих коэффициентов и с использованием структуры сетки в конструкторе оператора для всех ячеек рассчитываются разностные коэффициенты. В зависимости от того, с ячейками какого уровня (предыдущего, текущего или следующего) граничит данная ячейка, шаблон разностной схемы для неё может содержать от 7 до 25 точек.
Расчёты проведены на сетках 32 × 32 × 32 и 64 × 64 × 64 с 3, 4 и 5 уровнями адаптации. Решается нестационарное уравнение теплопроводности в кубе [0 ; 1]3 на промежутке времени t ∈ [0 ; 1]. Внутрь исходного куба помещается цилиндр вдоль оси z радиуса 0,2 с осью вдоль прямой (x=0,5; y=0,5). Внутри и вне цилиндра коэффициент диффузии κ полагается одинаковым и равным 1. Параметр a полагается постоянным и равным нулю. Возьмём точное решение в виде uexact(t;x,y,z) = e-t cos πx cos πy cos πz, тогда правая часть f(t; x,y,z) имеет вид f(t; x,y,z) = (3π2-1) ∙u(t; x,y,z).
Шаг по времени взят фиксированный. Тогда число явных шагов схемы ЛИ-М прямо пропорционально числу шагов по одному направлению (см. [14]), то есть при переходе с сетки 32 × 32 × 32 на сетку 64 × 64 × 64в стандартном варианте без адаптации время счёта возрастает в 16 раз. Ниже в таблицах приведено время счёта (в секундах) на указанных сетках при разном числе уровней адаптации.
Таблица 1. Время счёта при трёх уровнях адаптации. Основная сетка 32 × 32 × 32 64 × 64 × 64 «Старая» программа
4.99 90.31 gridmath
1.92 39.38 Таблица 2. Время счёта при четырёх уровнязх адаптации. Основная сетка 16 × 16 × 16 32 × 32 × 32 64 × 64 × 64 Последовательная версия программы примерно в 2,5 раза превосходит по производительности «старую» программу, написанную также на языке C++, но без использования данной библиотеки. Перекомпиляция для CUDA даёт ускорение ещё примерно в два раза. Также производилась декомпозиция сетки на несколько областей, каждая область обсчитывалась на отдельном узле, обмен граничными ячейками производился по MPI.
Не очень большое ускорение на CUDA объясняется спецификой задачи. CUDA даёт существенное ускорение, если удаётся сделать так, чтобы потоки с последовательными номерами обращались к соседним ячейкам памяти, что в случае локально-адаптивных сеток не выполняется. Этим же объясняется слабое ускорение при декомпозиции сеток. Основное время уходит на построение массива граничных ячеек, которые собираются по всей области. Также, в отличие от последовательной версии, появляется этап копирования данных из памяти CUDA в память host процессора. 5. Заключение
Написана операторная библиотека для решения задач математической физики на локальноадаптивных сетках, позволяющая компактно записывать в программах математические формулы, связывающие сеточные функции, определённые на таких сетках. Операторы в этой библиотеке превращены из абстрактных математических понятий в конкретные программные объекты. Перенос программ, использующих данную библиотеку, на графические ускорители CUDA, осуществляется простой перекомпиляцией.
С помощью библиотеки решена задача теплопроводности, заданная в виде параболического уравнения. Задача решается внутри прямоугольного параллелепипеда, в который помещено тело, на границе которого коэффициент диффузии может претерпевать сильный разрыв. Для повышения точности решения первоначально равномерная сетка измельчается в окрестности границы внутреннего тела. Перенос программы на CUDA (перекомпиляцией) даёт ускорение примерно в два раза по сравнению с последовательной версией.
Автор выражает ей искреннюю благодарность Феодоритовой О.Б. за консультации при переносе написанной ею программы, решающей параболическое уравнение на локальноадаптивных сетках, на данную библиотеку. Литература 1. The Blitz++ library. Official website. URL: http://blitz.sourceforge.net/ 2. NVidia CUDA.URL: http://www.nvidia.com/object/cuda_home_new.html 3. Intel Xeon Phi. URL:http://www.intel.com/content/www/us/en/processors/xeon/xeon-phidetail.html 4. DVM система. URL: http://www.keldysh.ru/dvm/ 5. OpenMP. URL: http://openmp.org 6. OpenACC URL: http://www.openacc-standard.org/ 8. TOP500 Supercomputer Sites. URL: http://www.top500.org 10. Гибридный вычислительный кластер K-100
URL: http://www.kiam.ru/MVS/resourses/k100.html Operator library for solving of mathematical physics problems on locally adaptive grids using CUDA Michael Krasnov This paper describes a library of operators to work on grid functions defined on three dimentional locally adaptive grids. The functions is designed so that the details of its implementation are hidden from the user, allowing effective implementation on machines with different implementation (hybrid, parallel, etc.), including CUDA. The library creation raises several goals: the approximation of the appearance of the program to the theoretical formulas, the relative ease of use, including the work on the graphics card with CUDA architecture, computational efficiency.