Разработка визуализатора движения жидкости

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

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

Введение

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

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

Данная тема вызывает интерес в силу следующих причин:

- сложность поведения моделируемого объекта, вызванная различными физическими явлениями, такими как конвекция, поверхностное натяжение и т. д., которые необходимо учитывать;

- жидкости могут взаимодействовать с объектами окружающей среды - стенками сосуда, куда помещена жидкость, скалами в море, плывущим кораблем;

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

- ограниченные вычислительные ресурсы, что особенно необходимо учитывать в симуляциях, которые должны работать в реальном времени.

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

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

Проблема вычислительных ресурсов, которые требуются для работы методов, основанных на вычислительной гидродинамике, может быть решена путем использования технологий GPGPU. Многие из существующих алгоритмов симуляции хорошо подвергаются распараллеливанию, а новые разрабатываются с учетом возможности применения этих технологий при реализации алгоритма. Таким образом, алгоритмы первой группы становятся хорошими кандидатами для применения в приложениях реального времени. В связи с этим в данной работе рассматривается алгоритм симуляции жидкости Position Based Fluids.

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

1. Изучить математические модели, применяемые для симуляции жидкости в реальном времени;

2. Выбрать математическую модель для симуляции жидкости;

3. Изучить алгоритм рендеринга жидкости в реальном времени;

4. Изучить методы использования технологииCUDA, а также необходимых вспомогательных библиотек;

5. Изучить методы использования библиотеки OpenGL версии не ниже 3.3 и необходимых вспомогательных библиотек;

6. Реализовать программу визуализации движения жидкости;

7. Разработать техническую документацию.

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

1. Описание математических моделей жидкости

Данная глава посвящена обзору существующих способов моделирования жидкости на основе методов вычислительной гидродинамики. Сперва приводится описание уравнений Навье-Стокса для несжимаемой жидкости. Затем производится обзор методов моделирования жидкости с точки зрения Эйлера и с точки зрения Лагранжа. Далее приведен обзор существующих на рынке аналогов программы. В конце приведены выводы по первой главе и обоснование выбора алгоритма Position Based Fluids для моделирования жидкости.

Уравнения Навье-Стокса

Уравнениями Навье-Стокса называется система из нелинейных дифференциальных уравнений в частных производных, описывающих движение жидкости. В частном случае, когда жидкость считается несжимаемой (плотность жидкости остается постоянной), в систему входят два уравнения: уравнение неразрывности (закон сохранения массы) и уравнение движения. Математически система записывается следующим образом [15]:

(1)

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

;

векторный дифференциальный оператор Лапласа, который определяется в трехмерном пространстве следующим образом:

;

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

Одним из способов моделирования жидкости является представление ее в виде поля, в котором каждому региону, содержащему жидкость, назначается набор свойств, таких как плотность, скорость, давление жидкости и т. д.[4]

трехмерный визуализация компьютерный

Рисунок 1 - Моделирование жидкости с точки зрения Эйлера

Этот подход также известен как сеточный метод (grid-based), поскольку на пространство симуляции накладывается сетка, разбивающая его на отдельные области, см. рисунок 1. Моделирование жидкости в целом производится путем рассмотрения процесса движения жидкости между ячейками. Дискретизация пространства может производиться различными способами. Сложная дискретизация пространства с использованием тетраэдрических и полигональных мешей обычно используется в программном обеспечении, предназначенном для решения инженерных задач, поскольку позволяет более точно обрабатывать поверхности со сложной геометрией [2]. Однако для целей компьютерной графики точность, обеспечиваемая подобным подходом, не является необходимой, поэтому в данной сфере отдается предпочтение более простым подходам.

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

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



Рисунок 2 - Моделирование жидкости с точки зрения Лагранжа

В отличие от сеточных методов, фиксирующих параметры жидкости в определенных точках пространства, которые не двигаются в течение симуляции, методы, основанные на подходе Лагранжа, привязывают параметры жидкости к частицам, которые могут перемещаться в пространстве, см. рисунок 2.Одним из таких методов является гидродинамика сглаженных частиц (SPH).

Изначально метод гидродинамики сглаженных частиц был разработан Гингольдом и Монаганом для осуществления астрофизических симуляций [10]. Разработчики хотели, чтобы их метод обеспечивал достаточную точность и был простым в применении. Впервые в компьютерной графике метод SPH был использован в работе Стама и Фиуме, посвященной визуализации огня и других газообразных явлений [17].Для моделирования жидкостей он был адаптирован Мюллером [11].

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

Функция ядра определяет, насколько большой вклад частица вносит в сумму из формулы (2).Данная функция должна быть:

· четной: ;

· нормализованной:

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

(3)

а для расчета лапласиана:

(4)

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

(5)

В правой части уравнения (5) находятся три поля плотности сил, соответствующих давлению, вязкости и внешним силам, которые совместно составляют поле плотности силы. Ускорение каждой частицы рассчитывается по следующей формуле:

(6)

При применении данного метода для оценки различных физических величин: плотности, давления и т. п. необходимо принять во внимание то, что формулы (2), (3), (4) могут не удовлетворять некоторым физическим принципам, например третьему закону Ньютона в случае оценки давления. Плотности силы давления для i-ой частицы оценивается по формуле, основанной на (4):

(7)

Силы давления, полученные по данной формуле, не будут симметричными, поскольку сила давления для частицы i будет рассчитана только на основе давления от частицы j (слагаемое с давлением от частицы i равно 0, т. к. градиент функции ядра, равен 0). Сила давления для частицы j также будет рассчитана только на основе давления от частицы i, а поскольку давления частиц iи j не обязательно одинаковы, силы, с которыми частицы воздействуют друг на друга не обязательно будут симметричными. Мюллер [11] устраняет эту проблему внося изменения в формулу (7):

(8)

Рисунок 3 - Пример работы метода SPH

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

Данный метод моделирования жидкости базируется на SPH. Основная идея метода - для каждой частицы осуществить предсказание плотности и на основе этого значения скорректировать давление так, что отклонение плотности каждой частицы от заданной плотности жидкости не будет превышать некоторого порогового значения [18].



Рисунок 4 - Алгоритм PCISPH

Алгоритм симуляции PCISPH приведен на рисунке 4. На каждой итерации цикла сперва осуществляется поиск соседей (строки 2-3), чтобы не осуществлять расчеты для частиц, которые не оказывают влияния. Затем вычисляются силы различной природы и скорости для всех частиц, при этом силы, возникающие в результате давления, полагаются равными 0(строки 4-7). Далее следует цикл предсказания и корректировки давления (строки 8-18), в котором для каждой частицы рассчитываются скорость, позиция, плотность и отклонение плотности от заданной для жидкости (строки 10-14). Полученное отклонение используется для коррекции давления по формуле:

(9)

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

 (10)

- функция ядра, j - радиус функции ядра.

Рисунок 5 - Пример работы PCISPH[18].

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

Данный метод основан фреймворке Position Based Dynamics (PBD), который широко используется в видеоиграх для симуляции физики твердых тел, ткани и т. п., поскольку обеспечивает достаточную скорость и визуальную правдоподобность симуляции [1]. Согласно этому методу, объект представляется в виде множества из N частиц и множества из M ограничений (constraints). Ограничением в PBDназывается уравнение или неравенство, ограничивающее движение тел относительно друг друга. Параметрами ограничений в PBD обычно выступают позиции и переменные, задающие ориентацию тела в пространстве. Кроме того, для ограничения может вводиться параметр жесткости, находящийся в пределах от 0 до 1, чтобы предоставить пользователю больший контроль над симуляцией.

Рисунок 6 - Алгоритм симуляции в PBD

Алгоритм симуляции в PBD [9]выглядит следующим образом, см. рисунок 6. В строках 1-3 производится инициализация всех частиц. Строки 4-18 представляют собой цикл симуляции. На каждой итерации предсказываются позиции частиц (строки 5-8) и генерируются ограничения (строки9-10). Затем производится разрешение ограничений (строки 11-13) в ходе которого предсказанные позиции частиц корректируются таким образом, что все ограничения выполняются. Затем производится обновление скоростей и позиций частиц (строки 14-17).

В алгоритме PBF в качестве ограничения используется следующее уравнение:



(11)

На его основе рассчитываются скорректированные позиции для всех частиц. Ограничение в формуле (11) призвано поддерживать плотность жидкости в симуляции постоянной и как следствие обеспечивать несжимаемость. Оценка различных физических величин и их градиентов в данном методе производится тем же образом, что и в методе SPH - при помощи функции ядра.

Рисунок 7 - Пример работы метода PBF

Метод PBF позволяет получить реалистичные визуальные результаты симуляции жидкости в реальном времени, см рисунок 7. Кроме того, он позволяет применять временные шаги большие, чем в методах основанных на SPH, при сохранении стабильности симуляции. Данный метод также решает проблему необходимости наличия определенного количества соседних частиц, свойственную методам SPH, вводя т. н. «искусственное давление». Это сохраняет визуальную правдоподобность моделируемой жидкости.

«Realflow» - программа, разработанная компанией «Next Limit Technologies», предназначена для создания моделей физических тел и симуляции различных физических явлений. Она может выступать как самостоятельное приложение, так и в роли плагина для 3D-редакторов, таких как Maya,Cinema 4D, 3dsMax и прочие. Изначально программа позволяла проводить только симуляции жидких тел, но на данный момент множество поддерживаемых физических тел расширилось и включает твердые, деформируемые тела, а также газы. Данная программа позволяет получить физически достоверные результаты при симуляции различных явлений ввиду применения физически корректных подходов. Однако это требует значительные объемы вычислений, поэтому симуляции производятся оффлайн, что не позволяет использовать «Realflow» в приложениях, работающих в реальном времени.

«Fluidsv.3» представляет собой симулятор жидкости с открытым исходным кодом. Для симуляции жидкости в нем применяются методы, основанные на SPH, в том числе описанные в данной главе:PCISPHи PBF.Применение технологий GPGPU таких как CUDA в данной программе позволяет производить симуляцию жидкости в реальном времени. Однако ввиду малого числа разработчиков, а также их нацеленности на разработку эффективных методов поиска соседей в алгоритмах, основанных на частицах, в программе отсутствует визуально правдоподобный рендеринг жидкости.

«FLIP Fluids» - платный плагин для 3D-редактора «Blender», позволяющий проводить симуляцию и осуществлять рендеринг жидкости. Для осуществления симуляций жидкости в данной программе применяются grid-based методы.«FLIPFluids» позволяет создавать такие эффекты, как водная пена, пузырьки воздуха в воде. Однако для достижения физической точности и получения необходимых визуальных результатов симуляция и рендеринг производятся офлайн, что не позволяет применять «FLIP Fluids» в приложениях реального времени.

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

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

2. Описание методов симуляции и визуализации жидкости

Position Based Fluids

В основе подходаPosition Based Fluids лежит фреймворк PBD и методы гидродинамики сглаженных частиц. Фреймворк PBD определяет структуру цикла симуляции и способ расчета позиций, а SPH - подход к оценке значений физических величин для частиц.

Позиции частиц рассчитываются на основе ограничений. Каждое ограничение состоит из:

· Размерности n;

· Функции ;

· Множества индексов;

· Опционального параметра жесткости;

· Тип ограничения: равенство или неравенство.

Считается, что ограничение, имеющее тип равенство выполняется, если , а имеющее тип неравенство, если . Параметр жесткости определяет, насколько сильно действует соответствующее ограничение.

Рисунок 8 - Общий алгоритм работы метода PBF

Общий алгоритм PBF представлен на рисунке 8. В строках 1-3 производится предсказание позиций частиц на основе скорости и временного шага. В строках 5-7 для каждой частицы производится поиск ее соседей, оказывающих влияние на данную частицу. Далее в цикле (строки 8-19) корректируются позиции частиц. После этого на основе новых позиций вычисляются скорости частиц и производится их постобработка (строки 21-22), а затем позиции частиц обновляются (строка 23).

Согласно PBD для каждой частицы вводится ограничение в виде следующего уравнения:

(12)

Ограничение , заданное формулой (12) моделирует несжимаемость жидкости и говорит о том, что плотность каждой частицы должна быть равна плотности жидкости.

Оценка плотности производится стандартным для методов SPH способом:

(13)

В качестве функции ядра для оценки плотности используется функция Poly6, график которой изображен на рисунке 9:

(14)

где h-радиус сглаживания функции ядра.



Рисунок 9 - График функции вдоль координатных осей при h = 1

Для оценки градиента используется функция Spiky, ее график показан на рисунке 10:

(15)

градиент которой имеет вид:

(16)



Рисунок 10 - График функций (непрерывный) и (пунктирный) вдоль координатных осей при h = 1.

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

(17)

Поиск осуществляется путем смещения вдоль градиента функции ограничения:

(18)

Градиент функции ограничения по отношению к частице с индексом рассчитывается по формуле:

 (19)

Формула (19) приобретает различный вид, в зависимости от того является частица с индексом соседом частицы с индексом или совпадает с -ой частицей:

(20)

Значение коэффициента может быть получено путем подстановки в уравнение (18):

(21)

Ввиду того, что является нелинейной, а также ее градиент на границе функции ядра исчезает, знаменатель в формуле (20) приводит к появлению нестабильности в симуляции, когда частицы находятся на границе функции ядра. Чтобы предотвратить это применяется техника constraint force mixing [16], согласно которой в формулу (18) добавляется член, призванный ослабить введенное ограничение:

(22)

Исходя из формулы (22) выражается следующим образом:

 (23)

Общее смещение для i-ой частицы может быть найдено как сумма смещений для i-ой частицы из всех ограничений, в которые входит i-ая частица:

(24)

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

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

Вместо этого в методе PBF применяется искусственное давление, которое задается при помощи функции ядра следующим образом:

(25)

Искусственное давление включается в формулу (24) следующим образом:

 (26)

Использование данного подхода позволяет поддерживать плотность для каждой частицы чуть меньшей, чем заданная плотность жидкости, а также создает эффект, похожий на поверхностное натяжение. Формула (26) является окончательной формулой для расчета смещения .

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

(27)

Затем рассчитывается корректирующая сила, на основе которой обновляются скорости:

(28)

В соответствии с методом XSPH скорости частиц изменяются следующим образом:

(29)

Применение данного метода позволяет добиться слаженного движения частиц жидкости.

Поиск соседних частиц

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

Uniformgridразбивает область, в которой происходит симуляция, на кубические ячейки одинакового размера. Длина стороны ячейки обычно полагается равной радиусу функции ядра, поэтому при симуляции в 3D для поиска всех соседей опрашиваются 27 ячеек[7].

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

1) Сounters- хранит число частиц, приписанных к соответствующей ячейке.

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

Затем каждая частица приписывается к соответствующей ячейке в отдельном потоке. Атомарный инкремент применяются чтобы корректно обновить счетчик частиц для ячейки, к которой будет отнесена частица. После этого индекс частицы записывается в соответствующий массив из CellParticles, см. рисунок 11.



Рисунок 11 - Пример построения uniformgridс помощью атомарных операций

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

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

Рисунок 12 - Пример построения uniformgridс помощью сортировки

Screen Space Fluid Rendering

Рендеринг в пространстве экрана (screen space rendering) это подход к рендерингу поверхности жидкости, для симуляции которой используются частицы [19]. Данный метод позволяет осуществлять рендеринг в реальном времени и регулировать соотношение между производительностью и качеством результата.



Рисунок 13 - Процесс рендеринга

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



Рисунок 14 - Пример карты глубины

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

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



Рисунок 15 - Пример карты толщины

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

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

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



Рисунок 16 - Пример карты сглаженной глубины

Для получения сглаженной карты глубины на рисунке 16 применяется двусторонний фильтр (bilateral filter). Данный фильтр выглядит следующим образом [14]:

(30)

Коэффициент нормализации рассчитывается как:

(31)

Двумерная функция ядра Гаусса вычисляется по следующей формуле:

(32)



Рисунок 17 - Пример карты нормалей

Пример карты нормалей приведен на рисунке 17. Цвет соответствует направлению нормали относительно наблюдателя: зеленый - вверх, синий -на наблюдателя, красный - вправо от наблюдателя.

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

(33)

Нормаль получается путем взятия векторного произведения между частными производными:

 (34)

.

Члены, содержащие полагаются равными 0, поскольку это упрощает вычисления и не влияет на результаты значительно.

Рисунок 18 - Пример финального изображения

Окончательный расчет цвета поверхности как на рисунке 18 производится по следующей формуле:

(35)

Уравнение Френеля позволяет определить доли отраженного и преломленного света. Оно используется, чтобы найти долю отраженного света, исходя из угла, под которым происходит наблюдение поверхности. Для осуществления рендеринга обычно применяется приближение данного уравнения, которое также называется приближением Шлика:

(36)

Базовая отражательная способность вычисляется по формуле:

(37)

Цвет преломленного света, прошедшего через жидкость, рассчитывается как:

(38)

Отражение и преломление могут быть реализованы с помощью техники под названием environmental mapping [8].

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

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

Рисунок 19 - Отражение луча от поверхности

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

Рисунок 20 - Применение техники environmental mapping'a для нахождения отраженного света

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

 (39)

На границе двух прозрачных сред свет преломляется согласно закону Снеллиуса:

(40)

Рисунок 21 - Применение техники environmental mapping'а для нахождения преломленного света

Вектор, обратный вектору падения преломленного света, см. рисунок 21, находится по формуле:

(41)

В данной главе было приведено детальное описание метода Position Based Fluids, выбранного для моделирования движения жидкости. Затем был рассмотрен способ поиска соседей, основанный на разбиении пространства на одинаковые кубические ячейки и описан процесс присвоения частиц к той или иной ячейке. В конце был разобран подход к визуализации жидкости под названиемScreen Space Fluid Rendering, описаны этапы построения карт визуальной глубины до поверхности и толщины жидкости, сглаживания карты визуальной глубины и извлечения из нее нормалей, а также построения финального изображения.

3. Разработка программы

Функциональные требования

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

Выбор ограниченного объема, в который помещается жидкость, перед запуском симуляции из предложенных вариантов;

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

Задание позиции жидкости внутриограниченного объема перед запуском симуляции;

Изменение параметров жидкости, таких как плотность, вязкость в течение симуляции;

Задание ускорения, действующего на жидкость вдоль координатных осей;

Задание числа итераций сглаживания поверхности жидкости;

Выбор цвета жидкости;

Задание коэффициента преломления жидкости;

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

Задание коэффициентов затухания света, проходящего сквозь жидкость.

Запуск симуляциив пошаговом и непрерывном режиме;

Постановка симуляции на паузу;

Перезапуск симуляции;

Перемещение камеры и изменение направления ее взгляда в 3D-сцене;

Вывод на экран 3D-сцены;

Вывод частоты кадров.

Технологии разработки

В качестве языка программирования для разработки данного приложения был выбран С++. Данный язык повсеместно используется для написания приложений, связанный с компьютерной графикой, о чем свидетельствует большое количество игровых движков и пакетов для работы с графикой, написанных на C++. Широкое применение этого языка в сфере компьютерной графики обосновано высокой скоростью исполнения программ, разработанных с его помощью, ввиду компиляции кода непосредственно в машинные инструкции. Также для C++ доступны различные интерфейсы программирования графических приложений: Direct3D, реализации спецификации OpenGL.Кроме того,C++ - один из языков, которыеCUDA позволяет использовать в качестве высокоуровневого языка программирования [3].

В качестве технологии GPGPU, с помощью которой реализована параллельная обработка частиц была выбрана CUDA. Одна из причин - уже упомянутая интеграция с языком C++. Помимо этого,CUDAобычно дает большую производительность на видеокартах NVIDIA по сравнению с OpenCL. Так же для CUDA существует множество библиотек, поддерживаемых NVIDIA и упрощающих разработку, что отличает CUDA от OpenGLComputeShaders. Примером библиотеки, которая используется в данной работе, является Thrust. Данная библиотека является аналогом библиотеки STL в C++ для CUDAи содержит реализации различных алгоритмов, приспособленных для параллельного выполнения на GPU.

Для осуществления рендеринга был выбран графический APIOpenGL, который, как и язык C++, широко применяется для разработки графических приложений. Вместе с OpenGL используется стандартный набор вспомогательных библиотек.Для инициализации окна, создания контекста OpenGL, обработки ввода и событий используется библиотека Graphics Library Frame Work (GLFW)[6]. Для работы с различными математическими объектами: векторами, матрицами, применяется библиотека OpenGL Mathematics (GLM)[13]. Для загрузки расширений и указателей на функции OpenGL во время исполнения программы используется библиотека Glad[5]. Для создания виджетов служит библиотека NanoGUI[12].

CUDA представляет собой платформу проведения GPGPU вычислений и модель программирования, созданные фирмой NVIDIA[3].В качестве платформы для вычислений CUDA предоставляет широкий спектр библиотек, содержащих готовые реализации различных алгоритмов, и расширения для нескольких языков программирования, основными из которых являются C и C++, см. рисунок 22.

Рисунок 22 - Предоставляемые библиотеки и поддерживаемые CUDA языки программирования

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

Рисунок 23 - Пример масштабирования программы для разного числа мультипроцессоров

Функции выполняемые на GPU в CUDA называются ядрами (kernels).Для того, чтобы определить функцию как ядро к ней необходимо добавить спецификатор __global__,

01

02

03

04

05

06

// Определение ядра в CUDA

__global__

voidCudaKernel(void* params)

{

…

}

Для вызова функции ядра применяется синтаксис конфигурации исполнения (executionconfiguration). С его помощью определяется размерность решетки (grid), в которую входят блоки потоков, и размерность самих блоков. Синтаксис выглядит как выражение <<<Dg, Db, Ns, S>>>, которое помещается между именем функции и списком ее параметров, описание переменных из данного выражения приведено в таблице 1.

Таблица 1 - Описание переменных для конфигурации исполнения ядра

Переменная	Тип	Назначение
Dg	dim3	Задает размерность решетки блоков.
Db	dim3	Задает размерность одного блока.
Ns	size_t	Задает размер разделяемой памяти в байтах, выделяемой для одного блока. По умолчанию 0.
S	cudaStream_t	Задает связанный поток. По умолчанию 0.

Тип данных dim3 представляет собой три беззнаковых целых числа. Если какие-то из трех чисел остаются незаданными, они полагаются равными 1. Возможность задавать размерность решетки и блока необходима для удобства работы с различными типами данных: векторами, матрицами, а также для проведения вычислений для некоторого объема. Примеры вызовов ядер с разной размерностью.

01

02

03

04

05

06

07

08

09

10

11

12// Вызовыядрав CUDA

constintblockSize = 128;

constintgridSize = 512;

//Размер решетки: (512, 1, 1), размер блока: (128, 1, 1)

CudaKernel<<<gridSize, blockSize>>>(params);

constdim3gridDim(256, 2);

constdim3blockDim(32, 4);

//Размеррешетки: (256, 2, 1), размерблока: (32, 4, 1)

CudaKernel<<<gridDim, blockDim>>>(params);

Пример иерархии потоков приведен на рисунке 24. В данном примере решетка имеет размерность (3, 2, 1), а каждый блок внутри решетки (4, 3, 1).

Рисунок 24 - Пример иерархии потоков[3].

Для функций, исполнение которых происходит на GPU, определены переменные, хранящие данные об индексах блоков и потоков, а также размерности блоков, потоков и решеток, в которых происходит исполнение этих функций, см таблицу 2.

Таблица 2 - Переменные с данными о размерностях и индексах

Переменная	Тип	Назначение
gridDim	dim3	Хранит размерность решетки.
blockDim	dim3	Хранит размерность блока.
blockIdx	uint3	Хранит индексы блока в решетке.
threadIdx	uint3	Хранит индексы потока в блоке.

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

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17 // Вычисление глобальных индексов потоков

// 1Dрешетка из 1Dблоков

__device__

int getGlobalIdx_1D_1D()

{

returnblockIdx.x * blockDim.x + threadIdx.x;

}

// 2Dрешетка из 2Dблоков

__device__

int getGlobalIdx_2D_2D()

{

intblockId = blockIdx.x + blockIdx.y * gridDim.x;

intthreadId = blockId * (blockDim.x * blockDim.y)

+ (threadIdx.y * blockDim.x) + threadIdx.x;

returnthreadId;

}

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



Рисунок 25 - Иерархия памяти в CUDA

Применение CUDA для реализации алгоритма Position Based Fluids необходимо для того, чтобы симуляция движения жидкости производилась в реальном времени.

Поскольку часть данных, которые будут подвергаться обработке при помощи CUDA, также будет использоваться для визуализации жидкости, возникает необходимость в обеспечении доступа CUDA к ним. К таким данным относятся позиции частиц, которые сохраняются в буферные объекты (bufferobject)OpenGL. Прежде чем буферный объект будет доступен для использования CUDA, он должен быть зарегистрирован в CUDA. Регистрация буферных объектов производится функцией cudaGraphicsGLRegisterBuffer, которая возвращает указатель на зарегистрированный объект cudaGraphicsResource*. Удаление ресурсов производится функцией cudaGraphicsGLUnregisterResource. Данные операции, обычно производятся для ресурса по одному разу.

01

02

03

04

05__host__?cudaError_t cudaGraphicsGLRegisterBuffer (

cudaGraphicsResource** resource, GLuint buffer, unsigned int flags)

__host__?cudaError_t cudaGraphicsUnregisterResource (

cudaGraphicsResource_t resource)

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

01

02

03

04

05__host__?cudaError_t cudaGraphicsMapResources(int count,

cudaGraphicsResource_t* resources, cudaStream_t stream = 0 )

__host__?cudaError_t cudaGraphicsUnmapResources(int count,

cudaGraphicsResource_t* resources, cudaStream_t stream = 0 )

Получить указатель на отображенный ресурс можно с помощью функции cudaGraphicsResourceGetMappedPointer, см. листинг 6.

01

02__host__?cudaError_t cudaGraphicsResourceGetMappedPointer(void** devPtr,

size_t* size, cudaGraphicsResource_t resource)

Далее приведено описание реализации алгоритма PositionBasedFluids с применением CUDA, см. рисунок 26.

Рисунок 26 - Соответствие функций этапам PBF

В фрагменте кода приведено ядро PredictPositions, в котором производится расчет новых позиций для каждой частицы.

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18__global__

voidPredictPositions(

constfloat3* positions,

float3* velocities,

float3* predictedPositions,

intparticleNumber,

float3gravityAcceleration,

floatdeltaTime)

{

int index = GetGlobalThreadIndex_1D_1D();

if (index >= particleNumber)

{

return;

}

velocities[index] +=gravityAcceleration*deltaTime;

predictedPositions[index] =positions[index] +velocities[index] *deltaTime;

}

Сперва производится расчет глобального индекса текущего потока (строка 10). Полученное значение используется для индексации частиц, поэтому дальше производится проверка на выход за границу массива частиц (строки 11-14). Это необходимо, поскольку число частиц не всегда равно числу потоков, осуществляющих их обработку. Потоки, для которых частиц нет, просто заканчивают свое выполнение. Далее обновляется скорость частицы на основе внешних сил, в роли которых здесь выступает сила тяжести (строка 16). Затем на основе обновленной скорости рассчитываются позиции, в которые должны попасть частицы (строка 17).

После предсказания позиций производится построение uniformgrid. С помощью функции transformиз библиотеки thrust, в которую передается функтор PositionToCellIdConverter рассчитываются хеш-значения ячеек для каждой частицы. Далее хеш-значения, позиции и скорости частиц сортируются с применением функции thrust::sort_by_key. После этого на основе отсортированного массива хешей производится поиск границ каждой ячейки в ядре FindCellStartEnd

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27voidFindCellStartEnd(unsignedint* cellIds, unsignedint* cellStarts,

unsignedint* cellEnds,intparticlesNumber)

{

extern__shared__unsignedintsharedCellIds[];

int index = GetGlobalThreadIndex_1D_1D();

if (index <particlesNumber)

{

sharedCellIds[threadIdx.x + 1] = cellIds[index];

if (index > 0 &&threadIdx.x == 0)

sharedCellIds[0] = cellIds[index - 1];

}

__syncthreads();

if (index <particlesNumber)

{

unsignedintcurrentCellId = sharedCellIds[threadIdx.x + 1];

unsignedintpreviousCellId = sharedCellIds[threadIdx.x];

if (index == 0 || currentCellId != previousCellId)

{

cellStarts[currentCellId] = index;

if (index != 0)

cellEnds[previousCellId] = index;

}

if (index == particlesNumber - 1)

cellEnds[currentCellId] = particlesNumber;

}

}

В строке 4 выделяется разделяемая память для блока потоков. Она служит для сохранения хеш-значений для частиц. Сперва каждый поток записывает в эту память хеш своей частицы (строка 9). Нулевой поток в блоке также записывает хеш предыдущей для него частицы из отсортированного массива (строки 10-11), если такая частица существует. Все потоки блока ожидают друг друга на барьере __syncthreads() - это гарантирует, что каждый из них записал хеш в разделяемую память. Далее каждый поток сравнивает хеш своей частицы с хешом предыдущей частицы, взятым из разделяемой памяти (строка 18). Затем, если поток обрабатывает нулевую частицу или хеши соседних частиц различны, индекс текущей частицы сохраняется как стартовый для ячейки текущей частицы (строка 20). При обработке ненулевой частицы ее индекс сохраняется как конечный для ячейки предыдущей частицы (строки 21-22). В строках 24-25 проверяется граничная ситуация, когда поток обрабатывает последнюю частицу.

Следующий этап - цикл коррекции позиций. В приведенном фрагменте кодаядра CalculateLambda из листинга 9производитсярасчетзначения из формулы (23) для всех частиц из одной ячейки с индексом cellId.

01

02

03

04

05

06

07

08

09

10

11

12

13for (int j = cellStarts[cellId]; j <cellEnds[cellId]; ++j)

{

float3positionDifference = positions[index] -positions[j];

floatsquaredPositionDifference = norm2(positionDifference);

densities[index] += poly6(squaredPositionDifference);

float3 gradient = spikyGradientKernel(positionDifference) *

restDensityInverse;

currentParticleGradient+=gradient;

if (index != j)

{

squaredGradientsSum += norm2(gradient);

}

}

Сперва вычисляется квадрат расстояние до соседней частицы и его квадрат (строки 3-4). Затем производится оценка плотности (строка 5) и градиента. После чего данный градиент добавляется к общему градиенту для текущей частицы (строка 8). Если же j-ая частица не совпадает с той, что обрабатывается потоком, квадрат длины градиента добавляется к переменной sqaredGradientsSum, которая соответствует .

В фрагменте кода на листинге 10 к sqaredGradientsSum добавляется квадрат длины градиента текущей частицы (строка 1).

Затем вычисляется значение функции ограничения по формуле (12) для текущей частицы (строка 2).Наконец вычисляется по формуле (23) (строка 3).

01

02

03 squaredGradientsSum += norm2(currentParticleGradient);

float constraint = densities[index] * restDensityInverse - 1.0f;

lambdas[index] = -constraint / (squaredGradientsSum + lambdaEpsilon);

Обновление позиций происходит в ядре CalculateNewPositions. Приведенный на листинге 11 фрагмент кода производит расчет смещения на основесоседних частиц из ячейки cellId согласно формуле (26). Полученное значение добавляется к общему смещению для текущей частицы (строка 6). После того как расчет deltaPosition будет завершен, производится проверка выхода частицы за границы симуляциии коррекция позиции при необходимости.

01

02

03

04

05

06

07

08for (int j = cellStarts[cellId]; j <cellEnds[cellId]; ++j)

{

float3positionDifference = positions[index] -positions[j];

floatlambdaCorr = correctionCoefficient *

powf(poly6Kernel(norm2(positionDifference)), n_corr);

deltaPosition += (lambdas[index] + lambdas[j] + lambdaCorr) *

spikyGradientKernel(positionDifference);

}

Далее производится расчет новых скоростей частиц при помощи функции thrust::transform и функтора VelocityUpdater. После этого производится постобработка скоростей.

В первую очередь находится векторзавихренностидля каждой частицы. Следующий фрагмент кода ядра CalculateVorticity производит вычисление завихренности для частиц-соседей из ячейки cellId на основе формулы (27) и прибавляет ее к переменной curl[index], которая соответствует вектору завихренности.

01

02

03

04

05

06

07 for (int j = cellStarts[cellId]; j <cellEnds[cellId]; ++j)

{

float3positionDifference = positions[index] -positions[j];

float3 gradient = spikyGradient(positionDifference);

float3velocityDifference = velocities[index] -velocities[j];

curl[index] +=cross(velocityDifference, gradient);

}

Далее вызывается функция ядра ApplyVorticityConfinement, корректирующая скорость, см. листинг 13. В строке 1 находится - градиент длины вектора завихренности, который затем нормализуется.После чего формула (28) применяется, чтобы вычислить (строка 9), используемую для коррекции скорости (строка 10).

01

02

03

04

05

06

07

08

09

10 float3eta = CalculateEta();

floatetaLength = length(eta);

float3normalizedEta{};

if (etaLength>epsilon)

{

normalizedEta= normalize(eta);

...