SPH法の実装(GPU実装含む)

SPH法をCPU,GPUで実装した例を示す．

近傍探索の効率化†

SPHなどの粒子法では近傍粒子からの影響を考慮して，支配方程式，物性値の近似を行うのだが，この近傍粒子探索に非常に計算時間がかかる．単純にすべての粒子の組み合わせについてそのユークリッド距離を計算した場合，粒子数Nの2乗のオーダで計算時間が増える．

近傍探索を効率化する手法として空間分割法がある．空間分割法とは探索したい物体がある空間を格子などで分割し，自身および隣接する分割領域に含まれる物体のみで距離を計算することで，計算時間を大幅に減らすことができる方法である．空間分割法の手順としては，

1. 物体の登録
   1. 空間を分割
   2. 各物体が所属する分割領域を計算
2. 近傍探索
   1. 探索中心座標が所属する領域を計算
   2. 自身と隣接する領域の物体をリストアップ
   3. リストアップされた物体との距離を算出

となり，大きく登録と探索の2つに分かれる．登録はシーン中の物体(粒子)位置が変化があった場合のみ行い，探索はその必要性があるときに逐次行われる．

空間分割法で問題となるのはどのように空間を分割するかである．分割方法の代表的な物としては，

• 等間隔格子(下図左)
• 八分木,四分木構造(下図右)
• kD木構造

などがあげられる．下の二つは分割領域を階層構造にして，領域に物体が含まれるかどうかで深さを変えられるため，等間隔格子に比べて探索効率が良い．ただ，物体の登録に時間がかかるため，同じ構造で多くの探索を行う場合，例えば，レイトレーシングなどでよく用いられている． 一方，粒子法では各ステップごとに粒子がダイナミックに動くため，毎ステップ領域への登録を行う必要がある． そのため，分割・登録が容易な等間隔格子がよく用いられる． 等間隔格子の場合は，階層構造を持つほかの領域分割法に比べて，登録作業を並列化しやすいという特徴も持つ．

　

GPUでの実装†

Particle Simulation using CUDA (PDF)(サンプルソースがNVIDIAのGPU Computing SDKに含まれる) を参考にしてSPH法へ拡張する．

最初にParticle Simulation using CUDAで用いられている近傍パーティクル探索手法について， 以下の図を例として説明する．

ステップの最初にパーティクルの登録作業を行う．登録手順は以下．

1. 各パーティクルのグリッドハッシュ値を計算．ハッシュ値は各グリッドセルでユニークで，かつ，セルの場所から計算できるものである(周囲グリッドを参照するため)．ここではグリッドセルの位置(i,j)から

   hash = i+j・nx

   とする．(nx,ny)はグリッドセルの数で，図2において各領域の左下に書いてある数字がハッシュ値である(nx=2)． ハッシュ値はGridParticleHash配列に格納する．また，後でソートしたインデックス値を格納する配列SortedIndexにパーティクルインデックス値を格納する．(GridParticleHash,SortedIndexのサイズ = パーティクル数)

   GridParticleHash	2	0	0	1	1	0	1	2
   SortedIndex	0	1	2	3	4	5	6	7

2. グリッドハッシュ値をキーとしてソートする．ソートには基数ソート(radix sort)が用いられている．CUDAを用いたradix sortについてはIPDPS '09の論文を参照．実装はCUDAサンプルにも含まれている．

   GridParticleHash	0	0	0	1	1	1	2	2
   SortedIndex	1	2	5	3	4	6	0	7

3. GridParticleHashで同じ値(同じ分割領域)が続く部分の始まりと終わりのインデックス(CellStartとCellEnd配列)を格納する． まず，配列(サイズ=グリッドセル数)を0xffffffffで初期化する．

   CellStart	0xff	0xff	0xff	0xff
   CellEnd	0xff	0xff	0xff	0xff

   ソート済みのGridParticleHashの各要素iに対して一つ前の要素(i-1)のハッシュ値を調べて，

   GridParticleHash[i] != GridParticleHash[i-1]

   ならば，iを始点としてCellStart[GridParticleHash[i]]に， i-1を一つ前のグリッドセルの終点としてCellEnd[GridParticleHash[i-1]]に格納する．

   CellStart	0	3	6	0xff
   CellEnd	2	5	7	0xff

   このときパーティクル位置などを格納した配列をSortedIndexを使ってソートする．

   SortedPos[i] = Pos[SortedIndex[i]];

近傍探索フェーズの手順は以下．

1. パーティクル位置 (or 任意の座標値) x のグリッドセルを算出
2. 周囲のグリッドを含めて以下の処理を実行
   1. グリッドハッシュ値hash計算
   2. CellStart,CellEndを参照して(CellStart[hash],CellEnd[hash])，グリッド内パーティクルの始まりと終わりのインデックス(start, end)を参照
   3. for(int j = start; j <= end; ++j) でループ
      1. 近傍候補パーティクル位置pjを参照

         pj=SortedPos[j]

      2. |pj-x|<hならそのパーティクルは近傍

SPHで用いる場合は，近傍探索フェースの最後で見つかった近傍パーティクルとの距離値を使ってカーネルを計算，積算していく． また，グリッド幅は有効半径hに基づいて決定する．

各パーティクルの密度値を計算する例を以下に示す．

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75

__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;
}
 

__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法では，陰関数で表現されたフィールドから等値面を抽出する手法である． 手順としては，

1. 立方体グリッド(Cube)を空間に配置
2. 各Cubeの辺上で陰関数値が0となる点を線型補間や二分法，ニュートン法などで算出
3. Cubeに含まれる点の配置と数から三角形メッシュを生成

グリッドに含まれるメッシュ数やメッシュの構成はテーブルを使って表す． より詳しい手順やソースコード，テーブルは以下のサイトを参照．

• http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/

GPUでの実装はCUDAでMC法(改造版)参照．

実装結果†

GeForceGTX580の環境でテストした． パーティクル数が約25000で，SPH計算だけなら200fps，メッシュ化も含めると80-90fpsぐらいになった． ちなみにCPUでも同様の方法で実装したところ，SPH計算だけで10fps，メッシュ化も含めると3fpsぐらいになった (メッシュ化には上記のGPU実装MC法を用いたが，GPUへのデータ転送時間がかなりかかっていると思われる)． 図3,4に実行結果を示す．

　　

ソースコード†

CUDA9.1用ソースコード†

Visual Studio 2015 + CUDA9.1の環境で動作確認したコードは以下．

rx_sph_v2.1.zip

• 必要ライブラリや説明は下記のCUDA5用のものとほぼ同じだが，デフォルトではプラットフォームがx64になっている点には注意(それにあわせてライブラリフォルダがshared/lib64，実行フォルダがbin64になっている)．
• 動作確認済み環境 : Visual Studio 2015 (Version 14.0.25431.01 Update3), CUDA9.1.85, freeglut3.0.0, GLEW2.0.0, boost1.62, ftgl2.1.3rc5(freetype2.7), FLTK1.3.3

CUDA5用ソースコード†

Visual Studio 2010 + CUDA5.0の環境で作成したコードは以下．

rx_sph_v2.zip

• ビルドするのに必要なライブラリ

  freeglut, GLEW, CUDA, boost, libpng, zlib, libjpeg, ftgl, FLTK

  各ライブラリについてはライブラリのインストールを参照．

• 簡単な説明＆注意事項
  • GUIツールキットとしてFLTKを用いている．
  • rx_fltk_glcanvas.hの59行目

    #define RX_USE_GPU

    をコメントアウトするとCPUで計算する．
  • Anisotropic Kernelによるメッシュ生成にも対応している．一応GPUでも実装しているが速度は遅め．
  • Visual Studioから実行する場合は，作業ディレクトリを"../../bin"にしておく．
  • CUDA5のSDKにあるヘッダファイル(helper_*)を用いているので，それらをshared/incにコピーするか，ヘッダがあるディレクトリをインクルードディレクトリとして設定すること．Windows 7だと標準では，

    C:\ProgramData\NVIDIA Corporation\CUDA Samples\v5.0\common\inc\

    にインストールされる．
  • "s"キー or Start/Stopボタンでシミュレーションスタート/ストップ．
  • F1〜F9キー or Sceneメニューからシーン切替可能．binフォルダに"sph_scene_*.cfg"という名前のファイルがあり， これがシーンを記述するファイル．デフォルトでは7つのシーンが入っている．
  • bin/cubeフォルダには何も入っていないが，屈折レンダリングしたい場合は，キューブマップ用の画像ファイル(6枚)を入れておくこと．ファイル名規則は*_negx.jpg, *_negy.jpg, *_negz.jpg, *_posx.jpg, *_posy.jpg, *_posz.jpg． サンプルとしてterragenで作成したものを以下に置く．
    cubemap_texture_example.zip
  • OBJ形式のファイルを読み込んで，固体オブジェクトにすることができる．テクスチャ付きでも大丈夫(png,jpgに対応)．
  • 研究用に作ったコードを削って軽くした物なので，いろいろ使っていない関数，変数，定義などが残っているかもしれない．