流體模擬：NeighborHood Search
流體模擬：NeighborHood Search
前言
領域搜尋（NeighborHood Search）對於流體模擬模擬是一個非常重要的問題，實現一個高效的領域搜尋可加速流體模擬，因此下面對我瞭解到的方法進行一定的介紹。
1. Uniform Grids（均一網格）
這個方法使用均勻的網格將模擬的空間細分為一個個大小均勻的單元組成的網格（對於粒子來說一般網格大小與粒子大小相等），這是最簡單的一種空間細分方法。我們在查詢領域粒子的時候，必須對相鄰單元中的粒子（即\(3\times 3 \times 3=27\)
使用這種方法的問題在於其中許多網格是空的，並無粒子存在，空間使用率低。
2. Spatial Hashing（空間雜湊）
採用空間雜湊的方法,首先我們根據粒子的位置計算對應的空間雜湊值，並將粒子對應的線性編號，將cell id與particle id儲存線上性陣列中，我們便可以得到一個cell id亂序的線性表；然後對cell id進行一個排序，這裡我們採用GPU的並行基數排序（可參考我的上一篇部落格：DirectX11:GPU基數排序 ），最後我們便可以得到一個cell id有序以及particle id按照cell id排序的兩個線性表。如下所示：
我們僅僅需要記錄每個空間單元cell記錄的粒子的起始與結束地址，然後根據當前網格得到其周圍27個空間單元的cell從而判斷其包含粒子與當前粒子是否接觸。
在這個方法中，所使用的雜湊函式如下：
\[\operatorname{hash}(\mathbf{c})=\left[\left(p_{1} i\right) \mathrm{XOR}\left(p_{2} j\right) \mathrm{XOR}\left(p_{3} k\right)\right] \bmod m \]其中\(p_1=73856093,p_2=19349663,p_3=83492791\)，輸入\(c=(i,j,k)\)
我所採用的方法也是基於空間雜湊的領域搜尋方法。
優點：稀疏表示，只儲存已填充的單元格
缺點：Cache命中不理想，是因為鄰居在記憶體中不緊湊；大質數（32位無符號整數可能不能準確表達）
3. HLSL核心程式碼
本人採用的是DirectX11中的computer shader進行實現，所有的緩衝區採用結構化緩衝區，假設你對以下的內容比較瞭解：
| DirectX11 With Windows SDK完整目錄 |
|---|
| 26 計算著色器：入門 |
| 深入理解與使用緩衝區資源 |
| 並行規約演算法 |
3.1計算雜湊值
//calculate hash
uint hashFunc(int3 key)
{
int p1 = 73856093 * key.x;
int p2 = 19349663 * key.y;
int p3 = 83492791 * key.z;
return p1 ^ p2 ^ p3;
}
uint GetCellHash(int3 cellPos)
{
uint hash = hashFunc(cellPos);
return hash % (2 * g_ParticleNums);
}
[numthreads(THREAD_NUM_X, 1, 1)]
void CS(uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x < g_ParticleNums)
{
float3 currPos = g_ParticlePosition[DTid.x];
int3 cellPos = floor((currPos - g_Bounds[0]) / g_CellSize);
g_cellHash[DTid.x] = GetCellHash(cellPos);
g_cellIndexs[DTid.x] = g_acticeIndexs[DTid.x];
}
else
{
g_cellHash[DTid.x] = -1;
g_cellIndexs[DTid.x] = -1;
}
}
g_Bounds為當前粒子位置各分量的最小/大值，可使用並行規約演算法中的上行階段（reduce階段）進行比較，從而得到最小的值。如下所示（求的是當前執行緒組資料的各分量的最小/大值）：
groupshared float3 boundmin[THREAD_BOUNDS_X];
groupshared float3 boundmax[THREAD_BOUNDS_X];
[numthreads(THREAD_BOUNDS_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID,uint GI:SV_GroupIndex,uint3 GId:SV_GroupID)
{
float3 currPos = float3(0.0f, 0.0f, 0.0f);
if (DTid.x < g_ParticleNums)
{
currPos = g_ParticlePosition[DTid.x];
}
else
{
currPos = g_ParticlePosition[g_ParticleNums - 1];
}
boundmin[GI] = currPos;
boundmax[GI] = currPos;
uint d = 0;
uint offset = 1;
uint totalNum = THREAD_BOUNDS_X;
//min
for (d = totalNum >> 1; d > 0; d >>= 1)
{
GroupMemoryBarrierWithGroupSync();
if (GI < d)
{
uint ai = offset * (2 * GI + 1) - 1;
uint bi = offset * (2 * GI + 2) - 1;
boundmin[bi] = min(boundmin[ai], boundmin[bi]);
boundmax[bi] = max(boundmax[ai], boundmax[bi]);
}
offset *= 2;
}
GroupMemoryBarrierWithGroupSync();
if (GI == 0)
{
float3 min = boundmin[totalNum - 1];
float3 max = boundmax[totalNum - 1];
g_ReadBoundsLower[GId.x] = min;
g_ReadBoundsUpper[GId.x] = max;
}
}
我們用粒子當前位置減去最小值，即可得到粒子在當前流體域中的網格位置資訊，從而計算其hash值並進行排序。
3.2 基數排序
詳情可見：DirectX11:GPU基數排序
3.3 當前網格的起始與結束地址
我們記錄每個空間單元cell記錄的粒子的起始與結束地址就變得非常簡單，只需將當前的cell id與排序好的線性表前一個的cell id進行比較，若不相同，則說明當前的粒子地址是當前粒子的所在空間單元cell的起始地址並且是前一個粒子所在單元空間cell的終止地址。
//load hashid at pos+1
groupshared uint sharedHash[THREAD_NUM_X+1];
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID,uint GI:SV_GroupIndex )
{
uint hashValue = 0;
if (DTid.x < g_ParticleNums)
{
hashValue = g_cellHash[DTid.x];
sharedHash[GI + 1] = hashValue;
if (GI == 0 && DTid.x > 0)
{
sharedHash[0] = g_cellHash[DTid.x - 1];
}
}
GroupMemoryBarrierWithGroupSync();
if (DTid.x < g_ParticleNums)
{
//If It's equal to the last one
if (DTid.x == 0 || hashValue != sharedHash[GI])
{
g_CellStart[hashValue] = DTid.x;
if (DTid.x > 0)
{
//the last cell end pos
g_CellEnd[sharedHash[GI]] = DTid.x;
}
}
//thread the last one output
if (DTid.x == g_ParticleNums - 1)
{
g_CellEnd[hashValue] = DTid.x + 1;
}
}
}
參考資料
[1] Green S. Particle simulation using cuda[J]. NVIDIA whitepaper, 2010, 6: 121-128.