![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
|
SPH法をCPU,GPUで実装した例を示す． 近傍探索の効率化†SPHなどの粒子法では近傍粒子からの影響を考慮して，支配方程式，物性値の近似を行うのだが，この近傍粒子探索に非常に計算時間がかかる．単純にすべての粒子の組み合わせについてそのユークリッド距離を計算した場合，粒子数Nの2乗のオーダで計算時間が増える． 近傍探索を効率化する手法として空間分割法がある．空間分割法とは探索したい物体がある空間を格子などで分割し，自身および隣接する分割領域に含まれる物体のみで距離を計算することで，計算時間を大幅に減らすことができる方法である．空間分割法の手順としては，
となり，大きく登録と探索の2つに分かれる．登録はシーン中の物体(粒子)位置が変化があった場合のみ行い，探索はその必要性があるときに逐次行われる． 空間分割法で問題となるのはどのように空間を分割するかである．分割方法の代表的な物としては，
などがあげられる．下の二つは分割領域を階層構造にして，領域に物体が含まれるかどうかで深さを変えられるため，等間隔格子に比べて探索効率が良い．ただ，物体の登録に時間がかかるため，同じ構造で多くの探索を行う場合，例えば，レイトレーシングなどでよく用いられている． 一方，粒子法では各ステップごとに粒子がダイナミックに動くため，毎ステップ領域への登録を行う必要がある． そのため，分割・登録が容易な等間隔格子がよく用いられる． 等間隔格子の場合は，階層構造を持つほかの領域分割法に比べて，登録作業を並列化しやすいという特徴も持つ．   図1 等間隔格子(左)と四分木(左)
GPUでの実装†Particle Simulation using CUDA (PDF)(サンプルソースがNVIDIAのGPU Computing SDKに含まれる) を参考にしてSPH法へ拡張する． 最初にParticle Simulation using CUDAで用いられている近傍パーティクル探索手法について， 以下の図を例として説明する．  図2 グリッドセルとパーティクル
ステップの最初にパーティクルの登録作業を行う．登録手順は以下．
近傍探索フェーズの手順は以下．
SPHで用いる場合は，近傍探索フェースの最後で見つかった近傍パーティクルとの距離値を使ってカーネルを計算，積算していく． また，グリッド幅は有効半径hに基づいて決定する． 各パーティクルの密度値を計算する例を以下に示す． /*!
* 与えられたセル内のパーティクルとの距離から密度を計算
* @param[in] gridPos グリッド位置
* @param[in] index パーティクルインデックス
* @param[in] pos 計算座標
* @param[in] cell パーティクルグリッドデータ
* @return セル内のパーティクルから計算した密度値
*/
__device__
float calDensityCell(int3 gridPos, uint i, float3 pos0, rxParticleCell cell)
{
uint gridHash = calcGridHash(gridPos);
// セル内のパーティクルのスタートインデックス
uint startIndex = cell.dCellStart[gridHash];
float h = params.EffectiveRadius;
float dens = 0.0f;
if(startIndex != 0xffffffff){ // セルが空でないかのチェック
// セル内のパーティクルで反復
uint endIndex = cell.dCellEnd[gridHash];
for(uint j = startIndex; j < endIndex; ++j){
//if(j == i) continue;
float3 pos1 = make_float3(cell.dSortedPos[j]);
float3 rij = pos0-pos1;
float r = length(rij);
if(r <= h){
float q = h*h-r*r;
dens += params.Mass*params.Wpoly6*q*q*q;
}
}
}
return dens;
}
/*!
* パーティクル密度計算(カーネル関数)
* @param[out] newDens パーティクル密度
* @param[out] newPres パーティクル圧力
* @param[in] cell パーティクルグリッドデータ
*/
__global__
void sphCalDensity(float* newDens, float* newPres, rxParticleCell cell)
{
// パーティクルインデックス
uint index = __mul24(blockIdx.x,blockDim.x)+threadIdx.x;
if(index >= cell.uNumParticles) return;
float3 pos = make_float3(cell.dSortedPos[index]); // パーティクル位置
float h = params.EffectiveRadius;
// パーティクル周囲のグリッド
int3 grid_pos0, grid_pos1;
grid_pos0 = calcGridPos(pos-make_float3(h));
grid_pos1 = calcGridPos(pos+make_float3(h));
// 周囲のグリッドも含めて近傍探索，密度計算
float dens = 0.0f;
for(int z = grid_pos0.z; z <= grid_pos1.z; ++z){
for(int y = grid_pos0.y; y <= grid_pos1.y; ++y){
for(int x = grid_pos0.x; x <= grid_pos1.x; ++x){
int3 n_grid_pos = make_int3(x, y, z);
dens += calDensityCell(n_grid_pos, index, pos, cell);
}
}
}
// 密度値を結果に書き込み
uint oIdx = cell.dSortedIndex[index];
newDens[oIdx] = dens;
}
デバイス関数calDensityCellはセル内のパーティクルを参照して，距離h以内ならカーネル関数を計算して密度値に積算していく． 表面メッシュ生成†パーティクルをそのままレンダリングしたのでは流体には見えない． 特に透明な液体をレンダリングする際には，液体表面の位置と法線を算出しなければならない． ここでは3DCGにおいて一般的に用いられている三角形メッシュを使って液体表面を表現する． 三角形メッシュによる液体表面抽出にはMC(Marching Cubes)法がよく用いられている． MC法では，陰関数で表現されたフィールドから等値面を抽出する手法である． 手順としては，
グリッドに含まれるメッシュ数やメッシュの構成はテーブルを使って表す． より詳しい手順やソースコード，テーブルは以下のサイトを参照． GPUでの実装はCUDAでMC法(改造版)参照． 実装結果†GeForceGTX580の環境でテストした． パーティクル数が約25000で，SPH計算だけなら200fps，メッシュ化も含めると80-90fpsぐらいになった． ちなみにCPUでも同様の方法で実装したところ，SPH計算だけで10fps，メッシュ化も含めると3fpsぐらいになった (メッシュ化には上記のGPU実装MC法を用いたが，GPUへのデータ転送時間がかなりかかっていると思われる)． 図3,4に実行結果を示す．    図3 左からパーティクル，表面メッシュ，屈折レンダリング
 図4 ポリゴンで表された谷形状内での流れ
ソースコード†CUDA9.1用ソースコード†Visual Studio 2015 + CUDA9.1の環境で動作確認したコードは以下．
CUDA5用ソースコード†Visual Studio 2010 + CUDA5.0の環境で作成したコードは以下．
|
cubemap_texture_example.zip 2214件
[詳細]
particle_grid.gif 3191件
[詳細]
particle_grid_2.gif 3291件
[詳細]
particle_octree.gif 3066件
[詳細]
rx_sph_v1.zip 1656件
[詳細]
rx_sph_v2.1.zip 2416件
[詳細]
rx_sph_v2.zip 3401件
[詳細]
sph_result1.jpg 4096件
[詳細]
sph_result2.jpg 3867件
[詳細]
sph_result3.jpg 3729件
[詳細]
sph_result4.jpg 3655件
[詳細]
sph_result4.png 1544件
[詳細]
