Моделирование распространения воздушных потоков и дыма с помощью графического процессора

Размещено на http://www.allbest.ru

Размещено на http://www.allbest.ru

Введение

На сегодняшний день в мире компьютерной сфере существует большое количество программного обеспечения, позволяющего создавать поражающие воображение спецэффекты для фильмов, мультфильмов, музыкальных клипов, новостных, рекламных и других видов видеороликов. Исторически все вычисления, связанные c моделирование спецэффектов, а также их визуализация и конечная обработка производилась на CPU.

Постепенно ситуация начала меняться и появились первые реализации фильтров пост-обработки, использующие GPU. Среди программного обеспечения, связанного с видео, такие фильтры одними из первых стали использовать графические ускорители. В силу естественности параллелизации GPU-реализации таких фильтров появились во всех профессиональных видеоредакторах, таких как, Adobe After Effects, Pinnacle Studio, Vegas.

Немного медленнее ситуация меняется среди программных систем, осуществляющих рендеринг, визуализацию готовых сцен. Первый коммерческий визуализатор, использующий GPU - NVIDIAOptiX - появился только в 2009 году, через 2 года после появления технологии CUDA, позволяющей проводить вычисления на GPU. Тем не менее на момент написания данного текста уже представлено несколько продуктов, способных осуществлять рендеринг, как полностью на GPU, так и в гибридном режиме - Chaos Group V-Ray RT, Octane, Arion, FurryballRT, Redshift, Cycles, moskito Render, Indigo Render, Thea Presto Render и т.д.

Еще медленнее ситуация изменяется среди программного обеспечения, позволяющего создавать спецэффекты, в т. ч. моделирующего физические явления: например, поведение жидкостей, газов и дыма. При этом на момент написания данного текста существующее ПО обладает рядом недостатков:

· Если сцена не помещается в память видеоускорителя, установленного на машине, то система моделирования, вообще, отказывается производить вычисления, а значит создавать спецэффекты для больших по размеру сцен становится чрезвычайно неудобно.

· Не существует доступных общественности приложений, позволяющих производить симуляцию спецэффектов, используя память нескольких GPU для одной задачи.

· Чаще всего GPU используются только для отображении примерного вида спецэффекта в окне предпросмотра, а не для финального рендеринга.

· Внешний вид эффектов, полученных при использовании GPU, отличается от вида эффектов, полученных с использованием CPU. В частности, финальная картинка, полученная на рендер-ферме с CPU слабо соответствует визуализации полученной в окне предпросмотра при использовании GPU.

При всем этом, несмотря на очевидные проблемы, художники по спецэффектам уже сейчас активно используют существующие плагины, например, ILM Plume, SideFX Houdini, Blender Cycles, поддерживающие видеоускорители, чтобы сэкономить время создания новой компьютерной Рисуноки. Еще большей пользы можно было бы добиться, если бы в широком доступе существовало ПО, способное производить не только предварительные расчеты, но и работать с большими сетками на GPU-кластерах, генерируя спецэффекты для финальной картинки. В моменты цейтнота и приближающегося срока сдачи того или иного фильма или ролика даже в крупных по размерах студиях художники подключают по сети свои персональные компьютеры к рендер-ферме, таким образом увеличивая суммарную производительность рендер-фермы и уменьшая временные затраты. При этом большое количество ресурсов в виде 3D-ускорителей простаивает вхолостую. Схожая ситуация наблюдается в маленьких ив домашних студиях. Там художник, обладая одним или несколькими мощными компьютерами лишь частично использует их ресурсы.

Существует множество причин сложившейся ситуации. Популярные комплексы по моделированию спецэффектов развиваются на протяжении многих лет, и портировать все наработки за короткое время с одной аппаратной платформы на другую невозможно (большая часть ПО, использующего GPU, написано с нуля и сильно отличается по возможностям конфигурации от CPU-аналогов). Средние по размеру студии не способны резко перейти на GPU-кластеры, не обладая гарантией возможности переноса большей части вычислений на видеокарты. Разработчики ПО не готовы производить продукт, который с большой вероятностью пока никто не купит. Возникает замкнутый круг.

На данный момент ПО способное эффективно задействовать вычислительные мощности GPU и опробованное на реальных проектах, существует лишь в нескольких особенно крупных студиях мира и недоступно широкой общественности, например, ILM Plume. Тем не менее, очевидно, что интерес к GPU-ориентированному ПО для создания спецэффектов растет, и с большой вероятностью рано или поздно настанет момент, когда видеоускорители повсеместно будут использоваться для создания спецэффектов.

1. Обзор существующих работ

Первые попытки смоделировать газ и дым в кинематографе были произведены в середине 80-х годов XX века. Одними из первых исследователей в этой области были Дж. Гарднер, который предложил математическую модель для моделирования и визуализации разных типовоблаков, а также У. Ривз, описавший создание спецэффектов с использованием частиц. Но работах того времени исследователи даже не пытались моделировать физическую суть эффектов.

Решение системы Навье-Стокса было слишком тяжелым для компьютеров того времени. Но дело было не только в произвольности - исследователи не знали, как управлять полученной субстанцией. Активно начали развиваться подходы, представляющие дым, как систему частиц с некоторой силой влияющих друг на друга. Такой подход обладает не слишком высоким качеством, но высокой производительностью, поэтому зачастую использовался в компьютерных играх. Подобные методы развиваются до сих пор. Например, предлагается алгоритм моделирования дыма, когда нужна высокая производительность, но качество не так важно.

Положение дел оставалось таким вплоть до появления статьи Н. Фостера, содержащей детальное описание алгоритма моделирования воды на основе решения уравнений Навье-Стокса. В дальнейшем при участии Дж. Стэма была опубликована серия работ, в которой при визуализации дыма среда также описывалась с помощью трёхмерных уравнений Навье-Стокса. Здесь жидкость и газ уже описывались, как эйлерова среда, но на некоторых этапах решения системы уравнений представлялись набором частичек. Поэтому такие алгоритмы стали называть полулагранжевыми. Алгоритм Стэма до сих пор считается каноническим, и на нем основано большинство последующих работ.

Для того, чтобы моделирование дыма было детальным нужно использовать особые процедуры, восстанавливающие завихрения, потерянные при расчётах грубыми методами, либо производить расчёты на очень больших сетках. Тут и встают вопросы памяти и производительности.

2. Постановка задачи

Моделируемая среда рассматривается как невязкий несжимаемый газ. Данные упрощения были сделаны из соображений, что расчёты будут проводиться на грубой сетке, и, таким образом, эффект вязкости в любом случае будет утерян. А так как скорость распространения дым, а заметно меньше скорости звука, то сжимаемостью также можно пренебречь. Рассматриваемый процесс описывается следующей системой трёхмерных уравнений Эйлера:

где и p - искомые вектор скорости и давление, а - известная внешняя сила, например, описанная далее сила плавучести.

Так как плотность считается константной и, в данном случае, равной единице, то в рассматриваемой системе имеется четыре неизвестных и четыре уравнения. На внешних границах имеющей форму параллепипеда области Gext, в которой будут проводиться расчёты, заданы естественные условия, позволяющие газу беспрепятственно входить и выходить из области. На границах объектов произвольной формы, расположенных внутри области Gint наложены условия не протекания, препятствующие попаданию дыма внутрь объектов. Начальное распределение искомых величин задаётся аниматором исходя из того, что за эффект он хочет воспроизвести, и в общем случае представляется заранее известными функциями U0 и p0, которые должны быть согласованы с граничными условиями.

Данная задача является основой для моделирования различных эффектов, таких как дым, огонь, облака - любых жидких и газообразных сред. Чтобы обеспечить моделирование именно дыма, в решаемую систему необходимо добавить ещё две неизвестные скалярные величины -- концентрацию c и температуру T дыма, который перемешивается и свободно распространяется вместе с основным газом. Именно эти величины в дальнейшем будут визуализироваться, в то время как газ считается бесцветным. Чтобы описать эту смесь из двух веществ, в введённую ранее систему добавляются два дополнительных уравнения адвекции, которые описывают процесс перемещения вещества в заданном поле скоростей и, как следствие, изменения его концентрации. После этого осталось дополнительно задать начальное распределение температуры и плотности дыма через также предоставленные аниматором функции:

И ввести обратную связь, чтобы тяжёлые клубы дыма опускались под действием силы тяжести, а потоки горячего воздуха восходили вверх - уже упомянутую силу Бузинеска. Для этого достаточно расписать функцию внешней силы следующим образом:

Где - сила тяжести, Tamb - температура окружающейсреды, а б и в - положительные константы, определяющие интенсивность процесса распространения дыма. Решив данную задачу для отрезка времени , мы получим данные, достаточные для наложения требуемого спецэффекта.

3. Метод решения

Расщепление.

Для решения введённой ранее системы уравнений предлагается задействовать вариацию метода расщепления, в рамках которого для перехода с i -ого на i + 1 временной слой должны быть выполнены три действия. Во-первых, исключив давление из уравнения сохранения момента импульса, необходимо вычислить предварительное значение вектора скорости :

Во-вторых, требуется обеспечить несжимаемость газа, для чего вычисляется актуальное значение давления путём решения уравнения Пуассона:

А в предварительный вектор скорости вносится соответствующая поправка:

В-третьих, необходимо найти новые значения для плотности и температуры дыма:

Адвекция.

Отдельно стоит остановиться на схеме численного решения уравнений адвекции, которые нужны для вычисления величин Если для этого использовать аппроксимацию конечной разностью, как это сделано в выкладках выше, то не удастся обеспечить абсолютную устойчивость метода, а наличие данного свойства является крайне важным в силу ограниченности вычислительных мощностей. Поэтому в рассматриваемом алгоритме используется полулагранжевый подход, в рамках которого на этапе решения уравнений адвекции -- и только на нём -- жидкость представляется набором частиц, расположенных в центрах ячеек специальной сетки (MACgrid, см Рисунок 1) и обладающих соответствующими характеристиками.

Например, в случае температуры это означает, что для узла с номером (I, j, l) требуется вычислить величину:

,

где - шаг по пространству. После этого значение на следующем временном слое будет определено как . Очевидно, что компоненты вектора практически всегда будут дробными, поэтому для вычисления искомой величины температуры необходимо провести интерполяцию значений из восьми ячеек сетки, расстояние которых до узла (i ?vx , j ?vy , l -vz) не превышает. Если же индексы выходят запределы расчётной области, то соответствующее значение определяется исходя из граничных условий.

Рисунок 1. Одна ячейка MAC-сетки

Важно заметить, что при переходе от этапа адвекции к этапу проекции мы будем считать, что все значение скоростей заданы на гранях, и направлению этих скоростей соответствуют орты, являющиеся нормалями граням к ячеек MAC решетки. Давление же для каждой частички нашей субстанции определенно в центрах ячеек.

Проекция.

Другим этапом решения исходной задачи является этап проекции. Он является наиболее ресурсоемким во всей предложенной схеме. Этот этап численно эквивалентен решению внешней задачи Неймана для уравнения Пуассона:

где:



Разностная аппроксимация данной задачи выглядит так:

А ее решение сводится к решению системы линейных алгебраических уравнений.

В данной работе не будет уделено внимание их реализации. Алгоритм работы этапа адвекции приведен лишь для введения понятия MAC-сетки, и описания требований к солверу для проекции. Основное внимание в этой работе будет уделено этапу проекции.

4. Выбор метода выполнения этапа проекции

Для того чтобы удостовериться в осмысленности использования графических ускорителей, и чтобы выбрать наилучший численный метод решения системы, было решено сравнить методы решения СЛАУ из нескольких готовых библиотек. Для GPU была выбрана библиотека AmgX, а для CPU - Hypre. Эти библиотеки обладают наиболее полным набором возможностей для соответствующих платформа. В частности, каждая из них имеет реализацию алгебраического многосеточного метода (AMG) и стабилизированного метода бисопряженных градиентов (BiCGStab) - наиболее производительных устойчивых современных вычислительных методов.

Тестирование было произведено на 3 типах сцен: без объектов внутри, с расположенной по центру сферой, со 125 равномерно распределенными по сцене сферами. Для валидации использовалась следующая задача математической физики:

Рисунок 2. Время достижения требуемой невязки: AMG, пустая область.

Рисунок 3. Время достижения требуемой невязки: AmgX область с выколотой сферой

На машине с Intel Xeon E5620, 32 GB RAM и NvidiaTesla K40c все тесты показали превосходство AmgX над Hypre. Ускорение на GPU достигало 20 раз. Кроме того, было выявлено превосходство AMG над BiCGStab на всех невязках для всех типов задач.

Полный многосеточный метод.

Как было показано наиболее эффективным методом решения уравнения на этапе проекции является многосеточный метод, вычисления в котором ведутся на сетках различной грубости. В AmgX нет возможности использовать информацию о геометрических свойствах сетки, но геометрическая вариация многосеточного метода эффективнее алгебраической, а также потребляет значительно меньше памяти. Следовательно, в свете вышеизложенной информации наиболее выгодным видится использование именно геометрического многосеточного метода.

Суть метода заключается в том, что уточнение приближенного решения на сетках более грубых по сравнению с целевой, позволяет увеличить скорость сходимости за счет уменьшение погрешности по всем частотам спектра разностного оператора.

Обозначим точное решение разностной задачи на сетке с шагом за , разностный оператор Лапласа за , погрешность решения на той же сетке за , а приближенно решение, за . Тогда . Подставим точное решение в исходную задачу.

Тогда:

где .

Полученное уравнение можно приближенно решить уже на более грубой сетке, а затем интерполировать его приближенное решение на изначальную сетку и сложить с приближенным решением исходной задачи, тем самым приблизив его к точному. Данную процедуру можно повторить для каждой сетки.

Имея это в виду, отметим, что при использовании некоторых численных методов на первых итерациях уменьшение вектора ошибки происходит независимо от величины шага сетки h. Такие методы называются сглаживателями.

Таким образом мы приходим к понятию V-цикла. Сначала выполняется несколько итераций сглаживателя на самой детализированной сетке. Затем вычисляется невязка, которая огрубляется на сетку следующая уровня и итерации сглаживания повторяются уже на ней. Итерации повторяются до тех пор, пока не будет достигнута самая грубая сетка. После этого результаты интерполируется на более точные сетки и складываются с приближенными решениями на них до тех пор, пока не будет достигнута самая точная сетка.

Рисунок 4.V-цикл

Способ, когда самая точная сетка не достигается с первого раза, а начиная с некоторого уровня вновь происходит огрубление, называется W-циклом. Последовательная комбинация V-циклов, которая начинается на самой грубой сетке и с каждым новым V-циклом достигает все более точных сеток, называется полным многосеточным методом (FullMultiGrid).

В данной работе используется именно полный многосеточный метод, так как он требует O(n) - в отличие от O(nlogn) для V- и W- циклов - операций для достижения требуемой точности, а значит асимптотически лучше, чем методы, состоящие из V- и W- циклов.

5. Реализация этапа проекции

Входные данные.

Как уже было сказано ранее, солвер для вычисления давления решает задачу численно эквивалентную решению уравнения Пуассона, однако внутри расчётной области могут быть произвольным образом разбросаны узлы, содержащие граничные условия: как Дирихле, так и Неймана.

Большая часть условий Дирихле равна нулю, т.к. обозначает те области пространства, в которых нет дыма. В том случае, если в поставленной задаче есть ненулевые условия Дирихле, они упаковываются в формат CSR и передаются в солвер.

С условиями Неймана все немного сложнее. У каждой ячейки MACgrid'а 6 нормалей, а значит, при подобном подходе нам придётся передавать в солвер 6 CSR-векторов с граничными условиями Неймана, а затем распаковывать их на GPU. В процессе вычисления так же придётся где-то хранить узлы «призраки», возникающие в результате расчёта условий Неймана.

Гораздо проще подставить формулы аппроксимации для «узлов-призраков» в изначальную систему разностных уравнений. В этом случае сами «узлы-призраки» будут вычеркнуты из системы, а значения граничных условий перейдут в .

Это можно выполнить следующим образом:

Предположим, что ячейка с индексом (i+1,j,k) заполнена твердым телом, тогда задано условие Неймана:

Подставим в исходное разностное уравнение. Получаем:

Заметим, что коэффициент при изменился, и он равен количеству невыколотых соседей . Запомним это - с коэффициентом мы будем работать на этапах сглаживания и подсчёта невязки. Так же заметим, что заданное значение условия Неймана вычитается из правой части. Это вычитание можно произвести уже сейчас на этапе подготовки входных массивов, что и выполняется.

Итого солвер получает на вход: маску, хранящую информацию о том, находится ли в данной точке дым, воздух или твердый объект, CSR сетку с ненулевыми значениями условий Дирихле и правую часть уравнения, исправленную в соответствии с заданными условиями Неймана.

Сглаживатель.

Первый этап работы многосеточного решателя - это сглаживание. В качестве сглаживателя разумно использовать схему «красное-черное» для метода Гаусса-Зейделя, позволяющую добиться эффективного распараллеливания на GPU.

Узлы с чётной суммой координат помещаются в одну «красную» группу, а узлы с нечётной суммой в другую «чёрную». Затем для каждой из групп поочередно выполняются итерации метода Якоби. Такое разбиение помогает CUDA-потокам при распараллеливании на GPU на каждом шаге либо только читать, либо только писать в узел, что делает вычисления детерминированными и значительно ускоряет производительность.

Таким образом, принимая во внимание главу о граничных условиях, нам необходимо выполнить:

Где верхние индексы обозначают итерацию, - значение правой части дополненное граничными условиями Неймана, а - количество соседних невыколотых узлов, причем , если (i,j,k) узел выколот.

Как уже было сказано, метод «красное-чёрное» позволяет добиться хорошего качества распараллеливания на GPU. Так как каждый из узлов одного цвета при выполнении итераций сглаживания зависит только от узлов другого цвета, внутри итерации для каждого из узлов выполняется либо чтение, либо запись, но никогда не обе операции. Данный подход позволяет достичь скорости метода Гаусса-Зейделя, использовав для расчетов на каждой итерации только одну сетку, а не две, как в случае с методом Якоби.



Рисунок 5. Красно-чёрная раскраска сетки. Срез на одном из слоев расчетной сетки. При расчёте черных узлов используются значения из красных и наоборот

Однако у данного подхода есть несовершенство при использовании на графических ускорителях - шаблон доступа к памяти. Графические ускорители устроены таким образом, что наибольшая производительность достигается, когда все потоки работают с участками памяти, расположенными подряд. По этой причине в данной реализации каждый расчётный слой хранится в линеаризованном массиве таком, что в нем сначала расположены только «красные» узлы, а затем только «чёрные». CUDA-ядра при этом работают на CUDA-grid'ах, вытянутых в линию. Таким образом на каждой итерации CUDA-ядра читают данные только из одного массива, а пишут в другой. Для того, чтобы не возникало путаницы, индекс в линеаризованном массиве рассчитывается специальной функцией по координатам узла, и такой подход используется для всех сеток внутри солвера, включая невязки, маски и правые части.



Рисунок 6. Линеаризация двумерного массива

Расчёт невязки

Расчёт невязки является операцией по реализации очень похожей на операцию сглаживания, поэтому работа с граничными условиями характерная для сглаживания характерна и для нее.

Невязка рассчитывается по формуле:

,

т.е.

Где количество соседей, считается так же, как в сглаживателе. Сокращение слагаемых происходит тем же образом.

Расчёт нормы невязки.

Чтобы уменьшить время расчёта, выполнение критерия останова проверялось не после V-цикла, а после выполнения итераций предварительного сглаживания на самой точной сетке. Такой подход, во-первых, позволяет сэкономить на расчёте невязки, так как в случае невыполнения критерия останова её все равно пришлось бы считать, а, во-вторых, на времени выполнения V-цикла. Так как операции предварительного сглаживания на точной сетке производятся гораздо быстрее, чем весь V-цикл, если критерий останова достигается за счет сглаживания до начала огрубления, то нет необходимости выполнять последующие вычисления.

В качестве метрики для проверки критерия останова использовалась норма Гёльдера:

Особенностью реализации подсчёта нормы на GPU является редукция переменных из warp'ов, которые, вообще говоря, могут не выполняться одновременно.

Был использован алгоритм, заключающийся в следующей последовательности действий:

1. Редукция внутри одного thread'а.

2. Редукция внутри одного warp'а.

3. Редукция внутри одного CUDA-блока.

4. Редукция внутри всего CUDA-grid'а.

Начиная с архитектуры Kepler, для видеокарт NVidia стало возможным получение в CUDA-потоке данных, лежащих в регистрах соседних CUDA-потоков.

Так как доступ к файлу регистров такой же быстрый, как доступ к L1-кэшу, данный подход в комбинации с методом сдваивания является наилучшим для редукции внутри одного warp'a.

На языке CUDA его можно описать так:

real warpReduceMax(realval) {

//Цикл по всем thread'ам warp'а, осуществляющий сдваивание

for (int offset = warpSize / 2; offset > 0; offset /= 2){

//Чтение данных из регистра другогоthread'а

real val2 = __shfl_down(val, offset);

val = (val> val2) ? val : val2;

}

returnval;

}

Данный подход не применим внутри всего CUDA-блока, т.к. чтение из соседних регистров возможно только внутри warp'a.

Для вычислений внутри блока в CUDA есть еще один тип памяти, работа с которой по производительности эквивалентна работе с L1 кэшем - shared (разделяемая память). Доступ к общей для нескольких thread'ов shared памяти возможен только внутри одного CUDA-блока, поэтому с ее использованием мы не можем произвести редукцию по всему grid'у, но можем выполнить ее внутри одного CUDA-block'a следующим образом:

1. Выделяем внутри для каждого блока массив в shared-памяти, вмещающий по одному числу с плавающей запятой, для каждого warp'a. Внутри одного блока может быть не более 1024 потоков - значит для того, чтобы вместить информацию из всех warp'ов достаточно 1024/32 =32 ячеек массива.

2. Нулевой поток каждого warp'а записывает результат внутренней редукции warp'а в отведенную для него ячейку.

3. Происходит синхронизация потоков.

4. Нулевой warp блока осуществляет редукцию указанного массива с помощью приведенного ранее алгоритма.

Теперь осталось произвести поблочную редукцию внутри grid'а. Здесь придётся воспользоваться глобальной памятью. В CUDA нет атомарной операции выбора максимального из двух чисел с плавающей запятой, но её можно реализовать с помощью атомарной операции сравнения и замены целых чисел, представив число с плавающей запятой, как целое:

#typedefunsignedlonglongint uint64

__device__double atomicMax(double * constaddress, constdoublevalue)

{

if (*address>= value)

{

return *address;

}

uint64 * const address_as_i = (uint64 *)address;

uint64 old = *address_as_i, assumed;

do

{

assumed = old;

if (__longlong_as_double(assumed) >= value)

{

break;

}

//выполняем сравнение и замену внутри атомарной операции

old = atomicCAS(address_as_i, assumed, __double_as_longlong(value));

} while (assumed != old);

return __longlong_as_double(assumed);

}

Операторы огрубления и пролонгации.

В данной работе используются полновзвешенный (full-weighted) 27-точечный оператор огрубления с билинейной интерполяцией (см Рисунок 7). Для пролонгации на более точные сетки так же используется оператор билинейной интерполяции (см Рисунок 8).

При этом, если узел, в который происходит пролонгация не является рабочим, т.е. является узлом заполненным воздухом или твердым телом, никакие данные в него не записываются, его значение остаётся равным нулю.



Рисунок 7. Веса узлов в операторе огрубления

Стоит отметить, что, если такой узел находится в зоне огрубления, то при замене значения этого узла на ноль. Оператор огрубления будет терять часть данных за счёт использования значения того узла, который не будет обновлен при пролонгации. Поэтому в данной работе оператор огрубления был модифицирован - знаменатели дробей, обозначающих вес того или иного узла, при огрублении рассчитываются динамически по следующему алгоритму:

1. Числители всех весов выставляются равными нулю. Общий знаменатель всех дробей выставляется равным нулю.

2. Если огрубляемый узел является рабочим, то его числитель приравнивается к единице, деленой на значение коэффициента пролонгации для данного узла, а общий делитель увеличивается на это значение.

3. Шаг 2 выполняется для всех точек огрубляемой области.

Такой оператор хорош тем, что не меняет сумму всех узлов при последовательном огрублении и пролонгации.

При пролонгации грубого решения на более точную сетку можно заметить, что сама по себе уточненная сетка не используется нигде, кроме сложения ее с приближенным решением на точной сетке. Объединения процессов пролонгации и сложения позволяет сэкономить память на хранении промежуточных данных, а также увеличить производительность за счет использования компилятором инструкций, производящей умножение и сложение за один такт.

Следует обратить внимание на то, что в каждый узел запись производится сразу несколькими потоками. Операция сложения в CUDA по умолчанию не является атомарной, а атомарная выполняется пропорционально количеству одновременных вызовов. Данную проблему можно обойти правильной последовательностью выполнения сложений - например, сначала сложить все узлы, находящиеся в верхних половинах пролонгируемых областей, затем в левых, затем произвести синхронизацию потоков и выполнить сложения для правых, а затем и нижних частей.

Кроме того, можно заметить, что часть узлов будет обновлено сглаживателем сразу же после пролонгации, а их значения не будут использованы. В данном случае это будут все узлы, помеченные для сглаживателя, как «красные» (см Рисунок 8), соответственно для увеличения производительности при пролонгации, вообще, не следует обновлять эти узлы.

Рисунок 8. Коэффициенты оператора пролонгации. Цветами отмечена "раскраска" узлов точной сетки.

Типы ячеек при огрублении определяются следующим образом:

1. Если хотя бы один узел в зоне огрубления содержит граничное условия Дирихле (ячейка заполнена воздухом либо задаёт константный поток газа на границе), то полученный узел считается так же имеющим условие Дирихле.

2. В противном случае, если в зоне огрубления есть хотя бы один рабочий узел (узел с газом), то полученный узел так же считается рабочим.

3. В остальных случаях узел считается выколотым (твердое тело).

Такой подход хорош тем, что на грубых сетках процент выколотых областей будет уменьшаться, а расчётных - увеличиваться (см. Рисунок 9). Это позволяет производить вычисления вплоть до самого грубого уровня.

Рисунок 9. Внешний вид маски расчётной области на сетках 129*129*129, 65*65*65 и 33*33*33 соответственно. Синим обозначены рабочие узлы, зеленым узлы с нулевыми граничными условиями Дирихле, красным - твердый объект с условиями Неймана на границе

сглаживание линеаризация двумерный массив

Взаимодействие между несколькими GPU.

До этого речь шла о расчётах только на одном графическом ускорителе, однако разработанное средство способно производить решение и на нескольких GPU.

Для этого используется технологияCUDAUnifiedVirtualAddressing. Начиная с версии 4.0 в CUDA появилась возможность создавать общее адресное пространство для всех вычислительных устройств в системе, то есть для того, чтобы скопировать данные с одного устройства на другое не нужно хранить информацию о том, что это за устройства, достаточно лишь знать указатели на источник и приемник данных, а функций cudaMemcpy сама понимает, как производить копирование. Более того, начиная с CUDA 6.0 и поколения Kepler, cudaMemcpy научилась при возможности автоматически использовать режим Peer-to-peer для прямого копирования данных с одного GPU на другой, избегая передачи данных через центральный процессор. Для активации этой возможности достаточно лишь включить в коде режим Peer-to-peer для каждого из используемых GPU. Именно этот механизм использован в данной работе.

Все данные равномерно распределяются между графическими ускорителями, при этом часть данных перекрывается двумя соседними GPU - у каждой пары соседних карт есть по два общих слоя данных (Рисунок 7).В сглаживателе каждый из ускорителей сначала производит вычисление тех значений, которые потребуются соседним ускорителям на следующих итерациях. Затем все невыколотые узлы упаковываются в буфер, размер которого равен максимально возможной длине пересылаемых данных. Так как маски общих узлов доступны каждому из соседних ускорителей, при обмене данным сосед легко может распаковать полученные данные.

Рисунок 10. Распределение данных между GPU. Синим обозначены данные на GPU0, зелёным - на GPU1, серым - выколотая область, а красным - область в которой постоянно происходит обмен данными между ускорителями

Выполняется асинхронный обмен данными с помощью функции cudaMemcpyAsync(). Пока происходит обмен данными, ускорители рассчитывают области, выделенные для них монопольно.

При использовании операторов огрубления и пролонгации, а также расчёта невязки вышеупомянутые буферы используются для хранения недостающих данных. Так, если одному из ускорителей требуется три слоя данных, а доступно только два, то недостающий слой записывается в буфер. Если доступен только один слой, то расчёты не производятся, так как их произведёт соседний ускоритель.

Результаты тестирования и валидации.

Валидация производилась с использованием модифицированной версии солвера, в которую кроме условий Дирихле так же передавались условия Неймана, что позволило численно решить задачу с существующим аналитическим решением.

В кубе

Аналитическое решение такой задачи находится с точностью до константы:

Однако если наложить дополнительное условие на внешние границы:

Получим аналитическое решение:

Рисунок 11. Зависимость погрешности решения от номера итерации. Сетка 129х129х129

Эту задачу удалось решить с помощью разработанного солвера. Зависимость погрешности решения от номера V-цикла отображена на Рисунке 11.

Также была проверена работоспособность программы на системе с двумя графическими ускорителями: NvidiaTeslaK20c и NvidiaTeslaK40c - и замерена производительность.

На двух картах было получено более чем полуторакратное ускорение решения.



Рисунок 12. Время решения задачи в зависимости от невязки на GPU. Сетка 257x257x257

Так же удалось полностью промоделировать решение системы Навье-Стокса, а затем произвести рендеринг итогового изображения с использованием BlenderCycles. Результат представлен на рисунке 13.

Рисунок 13. Рендер симуляции дыма

Заключение

Было разработано программное средство осуществляющее моделирование потоков дыма для систем компьютерной рисунков, позволяющее задействовать несколько графических ускорителей для наиболее ресурсоемкой части вычислений.

В частности, с помощью технологии NvidiaCUDA был реализован многосеточный решатель, осуществляющий этап проекции решения системы уравнений Навье-Стокса, в отличие от аналогов не ограниченный размером памяти лишь одного ускорителя, а эффективно задействующий все доступные системе устройства.

Полученный решатель был интегрирован в систему, реализующую полный цикл решения системы Навье-Стокса, позволяющую загружать данные о геометрии объектов, расположенных в моделируемой сцене, и выгружать данные о геометрии полученного дыма, что дало возможность встроить его в существующую систему рендеринга и получить итоговое изображение.

Литература

1. Lee V.W., Kim C., Chhugani J., Deisher M., Kim D., Nguyen A.D., Satish N., Smelyanskiy M., Chennupaty S., Hammarlund P., Singhal R. and Dubey P. Debunking the 100X GPU vs. CPU Myth: An Evaluation of Throughput Computing on CPU and GPU// Proceedings of the 37th annual international symposium on Computer architecture. Saint-Malo, France:2010.

2. Gardner G.Y. Visual Simulation of Clouds, Computer Graphics // SIGGRAPH 85 Conference Proceedings. 1985. P. 297-384.

3. Reeves W. Particle Systems - A Technique for Modeling a Class of Fuzzy Objects // ACM Transactions on Graphics.1983. 2. 2. P. 91-108.

4. Dong W., Zhang X., Zhang C. Smoke Simulation Based on Particle System in Virtual Environments //International Conference on Multimedia Communications. 2010.

5. Foster N., Metaxas D. Realistic Animation of Liquids, Graphical Models and Image Processing // Graphical Models and Image Processing archive. 1996. 58(5), P. 471-483.

6. Stam J. Stable Fluids // SIGGRAPH 99 Conference Proceedings, Annual Conference Series,1999. P. 121-128.

7. He S., Wong H., Pang W., Wong U. Real-time smoke simulation with improved turbulence by spatial adaptive vorticity confinement // Computer Animation & Virtual Worlds, Volume 22, Issue 2-3, 2011.P. 107-114.

8. Pfaff T., Thuerey N., Cohen J., Tariq S., Gross M. Scalable Fluid Simulation using Anisotropic Turbulence Particles // Proceedings of ACM SIGGRAPH Asia (Seoul, Korea, December 15-18, 2010), ACM Transactions on Graphics, vol. 29, no. 5, P. 174:1-174:8.

9. Fedkiw, R., Stam J., Jensen H. Visual Simulation of Smoke // In Proceedings of SIGGRAPH 2001, 2001.P. 15-22.

Размещено на Allbest.ru