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

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

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

Изменено: 24.02.2009

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

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

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





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

5. Поиск соседей

Оно же broad-phase (широкая фаза). Суть этого шага заключается в том, чтобы, получив на входе массив частиц, выдать на выходе такую структуру данных, которая потом позволила бы быстро получить доступ ко всем соседям определенной частицы. Для широкой фазы используется Spatial Hashing – разделение всего пространства на множество клеточек. Хэш каждой отдельно взятой частицы равен индексу клетки пространства, в которую она попадает. Для простоты я сделал количество клеточек строго фиксированным, а само пространство - ограниченным.

Тут надо сделать замечание по поводу хранения данных о частицах. Благодаря тому, что все частицы считаются одинаковыми, «личных» данных у них немного: позиция и скорость. Это векторные величины, для того, чтобы их запомнить нужно четыре скалярные переменные – по две на каждый ветор. То есть все нужные данные можно записать в один тексель текстуры: компоненты R и G содержат данные о позиции, а B и A – о скорости. Текстуру, содержащую данные обо всех частицах я буду называть «текстурой частиц» (напомню что, по сути, это просто массив частиц). Для этой текстуры был выбран 128-битный формат, т.к. остальные не обеспечивали необходимой точности.

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

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

Старая схема доступа к частицам по значению хэша
Рис. 7. Старая схема доступа к частицам по значению хэша

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

Новая схема
Рис. 8. Новая схема доступа к частицам по значению хэша

Теперь рассмотрим более подробно все этапы шага определения соседей. Для начала картинка:

Схема поиска соседей
Рис. 9. Схема поиска соседей

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

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

PS_OUTPUT ComputeHashPS ( float2 TextureUV : TEXCOORD0 )
{
	PS_OUTPUT Output;
	//позиция частицы будет находиться в part.xy
	float4 part = tex2D(ParticlesSampler, TextureUV);

	//Вычисление индекса текущей частицы
	//TexSize – размер текстур, хранящих индивидуальные данные частиц
	float id = TexSize*(floor(TexSize*TextureUV.y) + TextureUV.x);
	
	
	if (id < BodiesNumber)
		вычисляем хэш
		Output.RGBColor = float4(Hash(part.xy),id,0,0);
	else
		//хэш «лишних» частиц должен быть большим 
		//чтобы при сортировке они попали в конец таблицы
		Output.RGBColor = float4(CellSize*CellSize*4+1,70000,0,0);
	return Output;
}

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

Преобразование таблицы хэшей
Рис. 10. Преобразование таблицы хэшей

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

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

PS_OUTPUT AddressFindPS ( float2 TextureUV : TEXCOORD0 )
{
	PS_OUTPUT Output;
	float2 addrUV = floor(TextureUV*CellSize);
	int X = addrUV.y*CellSize + addrUV.x; //текущая ячейка

	float l = 0; 
	float u = TexSize*TexSize*part;
	bool found = false;
	float i = 0;

	float2 UV;
	int Ki;

	//обычный бинарный поиск
	while (u >= l && !found) 
	{
		i = (u + l)*0.5;
		UV = addrTrans(i);
		Ki = floor(tex2Dlod(HashSampler, float4(UV,0,0)).x);
		if (X < Ki){
			u = i-1;
		}
		else 
		{
			if(X > Ki){
				l = i+1;
			}
			else{
				found = true;
				break;//нужное значение найдено
			}
		}
	}

	float4 result = float4(0, 0, 0, 0);

	//если частица с таким хэшем обнаружена 
	if (found)
	{	
		//находим последнюю частицу с таким хэшем...
		int Ki = X;
		float f = i;
		while (X == Ki && f < TexSize*TexSize-1){
			f++;
			float2 UV = addrTrans(f);
			Ki = floor(tex2Dlod(HashSampler, float4(UV,0,0)).x);
		};

		//...и первую
		Ki = X;
		float b = i;
		while (X == Ki && b >= 0){
			b--;
			float2 UV = addrTrans(b);
			Ki = floor(tex2Dlod(HashSampler, float4(UV,0,0)).x);
		};
		//сохраняем результат
		result = float4((b+1), (f-b-1), 0, 0);
	};
	Output.RGBColor =  result;
	
	return Output;
}

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

PS_OUTPUT PermutateWithHashPS ( float2 TextureUV : TEXCOORD0 )
{
	PS_OUTPUT Output;

	float2 hash = tex2D(HashSampler, TextureUV);
	if (BodiesNumber > hash.y)
	{
		float2 UV = addrTrans(floor(hash.y));
		Output.RGBColor = tex2D(ParticlesSampler, UV);
	}	
	else
	{
	//лишние частицы ставим куда подальше :)
		Output.RGBColor = float4(wW*40,wH*40, 0, 0);
	}

	return Output;
}

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

6. Симуляция

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

9 ячеек
Рис. 11. 9 ячеек :)

Однако если вспомнить структуру адресной таблицы, то можно заметить, что данные, на которые ссылаются 1,2,3; 4,5,6 и 7,8,9 ячейки идут подряд.


Рис. 12. Цифры, указанные на ячейках в таблице частиц, указывают, каким ячейкам адресной таблицы они соответствуют

То есть, чтобы получить все необходимые частицы достаточно знать адреса 1, 4 и 7 ячеек и количество частиц в последующих 3-х ячейках. Поэтому я добавил еще один проход, который просто прибавляет к количеству частиц, лежащих в данной ячейке, частицы из двух следующих справа ячеек:

PS_OUTPUT ComposeAddr3PS ( float2 TextureUV : TEXCOORD0 )
{
	PS_OUTPUT Output;
	float4 current = tex2D(AddressSampler, TextureUV);
	float4 next1 = tex2D(AddressSampler, TextureUV + float2(invCellSize,0));
	float4 next2 = tex2D(AddressSampler, TextureUV + float2(invCellSize*2,0));
	Output.RGBColor = float4(current.x,current.y + next1.y + next2.y,0,0);

	if (current.x == 0.0f)
	{
		if (next1.x == 0.0f)
			Output.RGBColor.x = next2.x;
		else
			Output.RGBColor.x = next1.x;
	}	
	return Output;
}

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

float4 parti = tex2D(ParticlesSampler, OriginalUV); //данные о текущей частице
float2 pos = parti.xy;
	
float step = invCellSize;
float2 AaddrUV = pos*invwWH + half2;
float2 addrUV = (floor(AaddrUV*CellSize) + half2)*step;

addr[0] = tex2D(AddressSampler, addrUV + float2(-step,0));
addr[1] = tex2D(AddressSampler, addrUV + float2(-step,step));
addr[2] = tex2D(AddressSampler, addrUV + float2(-step,-step));
		
for (int j = 0; j < 3; j++)
{
	int i = 0;
	float2 partUV = addrTrans(addr[j].x); //координаты «соседа»
	while (i < addr[j].y)
	{
	//тут тело цикла

		partUV.x += invTexSize;
		if (partUV.x >= 1.0f)
			partUV += float2(-1.0f,invTexSize);
		i++;
	};
}; 

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

Для вычисления плотности достаточно просто перенести в тело цикла известные формулы SPH:

float2 partj = tex2Dlod(ParticlesSampler, float4(partUV,0,0)).xy;
float2 dist = (pos - partj);
float dsq = dot(dist,dist);
			
if (mR2 > dsq)
{
	float c = (mR2 - dsq);
	sum += c*c*c;			
} 

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

Ro = sum*sph_pmass*DensKern*d6;
P = (Ro - sph_restdensity)*sph_intstiff;
//добавление константы позволяет не бояться переполнения при нулевой плотности
Ro = 1.0f/(Ro+0.0000001); 

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

float2 tmpUV = floor((OriginalUV)*TexSize);
int X = (tmpUV.y*TexSize + tmpUV.x); //индекс текушей частицы
int dif = X - BodiesNumber;

if (dif >= 0 && dif < AddParticles) //если нужно добавить частиц
{
	float offset = -(dif - AddParticles*0.5)*0.15;
	float2 newpos = Emitter.xy + offset*float2(Emitter.w,-Emitter.z); //координаты

	float fi = X*0.5*5.0/AddParticles;
	float2 newvel = float2(4*sin(fi)+20*sin(fi/40.0)*EmitterMagn,25.0); // скорость
	float2x2 Rotate = {Emitter.w, Emitter.z, -Emitter.z, Emitter.w};
	newvel = mul(Rotate,newvel); //скорость с учетом поворота
		
	Output.RGBColor = float4(newpos,newvel);
	return Output;
}

Далее идет расчет сил взаимодействия частиц. Тут все очень похоже на предыдущий шаг, просто применяются формулы SPH. Умножение на все константные составляющие производится вне цикла, к целях оптимизации. Кстати говоря, именно цикл по всем соседям является узкими местом алгоритма. Переменная float2x2ce содержит данные о силах, вызванных давлениям в компонентах xy и данные о силе вязкости в zw.

float4 partj = tex2Dlod(ParticlesSampler, float4(partUV,0,0));
float4 dist = (parti - partj); // dist.xy – расстояние, dist.zw – разность скоростей
float r2 = dot(dist.xy,dist.xy) + 0.0001f;
			
if ( mR2 > r2)
{
float2 densj = tex2Dlod(DensTextureSampler, float4(partUV,0,0)).xy;
float r = sqrt(r2);
float c = (r - mR);
float pterm = c * (dens.y + densj.y)/(r);
float2 dfloat2x2ce = float2(pterm,densj.x); 
force += dist*dforce.xxyy*c; //обе силы обрабатываются одновременно
//тот же код, но с разделением по роду сил 
//force.xy += dist.xy*c*pterm;
//force.zw += dist.zw*c*densj.x; 
};

...

force.xy *= -0.5f*PressKern*0.0007f*d*d; //постоянные множители для давления
force.zw *= VescKern*sph_visc*d*d; //...и для вязкости

Следующий шаг – взаимодействие со статической геометрией. Я использовал наиболее простой способ: Статическая геометрия рендерится в карту нормалей (в моем случае - вручную), полученная карта используется для приложения сил к частицам. Сила, с которой частица выталкивается вдоль нормали, зависит от ее скорости и давления. Зависимость силы от глубины проникновения не заложена явно, но можно добиться похожего эффекта, если сделать нормали на краю геометрии более «вертикальными», а по мере углубления – более наклонными. Из x- и y-координатам нормали, полученной с помощью текстурной выборки нужно вычесть небольшую поправку, равную 0.5/255.0, в противном случае, сила выталкивания начнет действовать на те частицы, которые на самом деле не касаются геометрии.

if (bound)
{
	float2 boundUV = (floor(AaddrUV*bound_size) + half2)/bound_size;
	float2 bNormal = tex2D(BoundSampler, boundUV) - half2*(1.0 + 1.0/255);
	bNormal.y *= -1;
		
	float2 v = float2(0,0);
	if (vel.x*bNormal.x < 0) v.x = abs(vel.x);
	if (vel.y*bNormal.y < 0) v.y = abs(vel.y);
		
	float mult = 0.003f/dens.x;
	//сила пропорциональна величине xy компонент нормали
	force.xy += bNormal*bound_k*mult*mult*0.8 + v*bNormal*bound_kv*mult*3;
}

Следующим обрабатывается аттрактор – притяжение или отталкивание частиц с помощью мыши. Координаты мыши на экране транслируются в координаты мыши в мире и передаются в шейдер:

D3DXVECTOR3 v_pos(mousePoint.x*1.0f/BBw - 0.5,0.5 - mousePoint.y*1.0f/BBh,0.0);
D3DXVECTOR4 v_out;
D3DXVec3Transform(&v_out,&v_pos,&mWorldViewProjectionInv);
v_out *= v_out.w;
V( g_pPhysics->SetfloatArray(pHanldeArr[h_p_cur_position],v_out,2));

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

if (attractor)
{
	float2 attr_dist = cur_position - pos;
	float c = 1.0/(length(attr_dist)+0.1f);
	force.xy -= attr_dist*d*attractor_force*c*c*c*attractor;
}

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

float2 accel = (force.xy + force.zw)*dens.x*sph_pmass*10;
accel += SS_accel*0.5 + float2(0,g*gm);

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

float max_acc = 50.0f;
if (dot(accel,accel) > max_acc*max_acc)
	accel = normalize(accel)*max_acc;
	
FloorConstraintEX(pos, vel, accel);

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

void FloorConstraintEX(float2 pos, float2 vel, inout float2 accel)
{
	float diff = 0;
	float d = sph_simscale;
	float r2 = 2 * sph_pradius;

	float4 diff4 = r2 + (float4(-pos.y,pos.y,-pos.x,pos.x) - float4(wH,wH,wW,wW)*0.5)*d;
	diff4 = max(diff4,0);
	float4 extdamp4 = mul(fNormals,vel*d);

	float4 adj = diff4*sph_extstiff - sph_extdamp*extdamp4*step(0.0001,diff4);
	accel += mul(adj,fNormals)*0.01f;
}

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

vel += accel*dt*10/d;
pos += vel*dt;
FloorConstraint(pos, vel);
		
Output.RGBColor = float4(pos,vel);
return Output;

Все, дело сделано.

7. Диффузия

Чтобы не мешать все в одну кучу я решил выделить расчет диффузии в отдельный шаг.

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

Координаты данных о частице и о цвете частицы совпадают
Рис. 13. Координаты данных о частице и о цвете частицы совпадают

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

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

float2 tmpUV = floor((OriginalUV)*TexSize);
int X = (tmpUV.y*TexSize + tmpUV.x);
int dif = X - BodiesNumber;

if (dif >= 0 && dif < AddParticles)
{
	Output.RGBColor = pDiffuse;
	return Output;
}

Далее опять пробегаем циклом по всем соседям, как уже два раза делали ранее. В цикле же смотрим расстояние и цвет соседа. В зависимости от расстояния берем «частичку» его цвета и прибавляем к цвету нашей частицы. И, чтобы не нарушить баланс, вычитаем из нее такое же количество ее собственного цвета. Где-то в соседнем потоке нужной частице прибавят то, что мы тут вычитаем, так что все нормально.

float4 parti = tex2D(ParticlesSampler, OriginalUV);
float4 color = tex2D(ColormapSampler, OriginalUV);
float4 color1 = color; //сохраняем исходный цвет

...

//тело цикла
float4 partj = tex2Dlod(ParticlesSampler, float4(partUV,0,0));
float4 dist = (parti - partj);
float r2 = dot(dist.xy,dist.xy);
		
if ( mR2 > r2)
{
float4 colorj = tex2Dlod(ColormapSampler, float4(partUV,0,0));
float r = sqrt(r2);
	float c = (1 - r/mR);
	color -= color1*c*DiffK;
	color += colorj*c*DiffK;
};

...

Output.RGBColor = color;
return Output;

Единственное, что тут нужно иметь в виду, если использовать 32-битный цвет для частиц, то диффузия достаточно быстро прекратится, так как небольшие изменения цвета будут проигнорированы из-за значительной погрешности. Поэтому я использовал 64-битный формат.

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

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