Бесплатная реклама

Автор: Лев Симонов «LEVel»

Опубликовано: 24.02.2009

Изменено: 24.02.2009

Постоянная ссылка

Комментарии [4]

Моделирование жидкости в 2D с использованием GPU


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



Страницы: 1 2 3

Содержание

1. Пространное вступление
2. Общее описание техники
3. Smoothed Particle Hydrodynamics
4. Пара слов о GPGPU
5. Поиск соседей
6. Симуляция
7. Диффузия
8. Рендеринг
9. Демонстрация (исходный код)
10. Приложение. Сортировка методом bitonic merge sort на GPU

GPU Water

1. Пространное вступление

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

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

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

Почему 2D? Потому что с точки зрения математической модели 2d и 3d – практически одно и то же, только у векторов другая размерность. А вот рендерить трехмерную воду – задача не из легких, не говоря уже о том, чтобы сделать результат рендеринга красивым и реалистичным. В 2d же все необходимое делается за 1-2 вечера с раскачкой и перекурами. Тема рисования результата вычислений, кстати, тоже будет отдельно немного освещена.

Почему GPU? Потому что, ну… это круто! Нет, на самом деле, я не знаю, насколько оправдан перенос вычислений на видеокарту в данной конкретной задаче. Для того чтобы узнать это, нужно было бы сделать вдвое больше работы, осуществив все то же самое на CPU. В любом случае, моей целью было научиться считать на GPU, что и было, в той или иной степени, удачно осуществлено. Приобретенные знания, наработки и опыт окажутся полезными при проведении вычислений на графическом чипе в других задачах, где это может дать большой выигрыш в производительности. Я же постараюсь поделиться полученным опытом с вами.

Почему не CUDA? Все вычисления производятся в духе «старой школы» GPGPU, то есть непосредственно с помощью графического API (в моем случае - DirectX). Это оправданный шаг, потому что использование CUDA автоматически выводит из игры все видеокарты, кроме nvidia 8ххх, 9ххх и 2хх серий. Мой же вариант работает почти везде, где есть шейдеры версии 3.0. Если, скажем, разрабатывается компьютерная игра, а не научная работа, то пренебрегать поддержкой такой большей доли графических адаптеров просто нельзя.

Теперь, когда все встало на свои места, можно приступить к делу.

2. Общее описание техники

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

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

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

Общая схема
Рис. 1. Общая схема

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

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

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

Наконец, рисуется интерфейс и справочная информация, эта задача лежит на плечах DXUT.

3. Smoothed Particle Hydrodynamics

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

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

где A – значение некоторой скалярной величины, m – масса частицы, – плотность.

Градиент и лапласиан A берутся просто как

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

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

Нам нужно вычислить три составляющие плотности силы: давление, вязкость и внешние силы (в данном уравнении - гравитация). Для вычисления плотности силы давления используется следующая формула:

здесь p – давление, вычисленное по формуле

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

Плотность силы вязкости вычисляется по формуле:

Конечное ускорение частицы можно представить в виде:

где - сумма внешних сил, действующих на частицу.

В своей программе я использовал такие же коэффициенты и функции сглаживания, как и в работе Fluids v.1.

4. Пара слов о GPGPU

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

Основное, на мой взгляд, отличие вычислений на GPU от CPU заключается в работе с памятью. Во-первых, нельзя писать данные в произвольную ячейку, во-вторых, количество обращений к памяти необходимо минимизировать, так как пропускная способность ограничена и может легко стать узким местом. Исходные данные передаются в GPU в виде текстур, играющих роль массивов, и констант. А результаты вычислений оказываются в рендертаргетах, опять же, играющих роль массивов. Выглядит это примерно так:

Вычисления на GPU
Рис. 2. Вычисления на GPU

Видно, что для каждого текселя рендертаргета происходит выполнение кода пиксельного шейдера. На языке CPU это выглядит как цикл по всем текселям с выполнением определенного кода.

Кстати о массивах. Обычно данные удобно хранить в одномерных массивах, но размерность одномерных текстур в современных видеокартах ограничена, и, как правило, не может превышать 4096 элементов. В то же время, двумерные текстуры могут достигать размера 4096х4096. Очевидно, что одномерные массивы лучше запаковать в двумерные текстуры таким образом:

Трансляция одномерного массива в двумерный
Рис. 3. Трансляция одномерного массива в двумерный

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

float2 addrTrans(float address1D)
{ 
	float2 res;
	res.x = modf(address1D*invTexSize, res.y);
	res.y *= invTexSize;
	return res;
}

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

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

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

Кстати о квадах. Для простоты я решил не менять вершины динамически, а просто использовал один универсальный квад:

Квад
Рис. 5. Квад. Черные числа – координаты вершин, зеленые – текстурные координаты

Если бы нужно было, скажем, просто скопировать одну текстуру в другую, то код шейдеров был бы таким:

void SAQuadRenderVS( float4 Pos : POSITION,
			float2 UV : TEXCOORD,
                    out float4 oPos : POSITION,
                    out float2 oUV : TEXCOORD0)
{
    oPos = Pos;
    oUV = UV;
}

PS_OUTPUT ComputeHashPS ( float2 TextureUV : TEXCOORD0 )
{
	PS_OUTPUT Output;
	Output.RGBColor = tex2D(Sampler, TextureUV);
	return Output;
}

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

Задача обработки произвольной области данных передо мной так и не встала, но в целях оптимизации я добавил функцию «стягивания» квада к верхней границе текстуры, что соответствует выполнению цикла не по всем элементам массива, а только по первым N*n:

«Стягивание» рабочей области с разными коэффициентами
Рис. 6. «Стягивание» рабочей области с разными коэффициентами

Делается такое стягивание простым смещением координат в вершинном шейдере:

void SAQuadRenderPartVS( float4 Pos : POSITION,
			 float2 UV : TEXCOORD,
                     out float4 oPos : POSITION,
                     out float2 oUV : TEXCOORD0)
{
Pos.y = Pos.y*n + 1.0f - n;
UV.y *= n;

oPos = Pos;
oUV = UV;
}

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

Страницы: 1 2 3