Curl Noise追記

GPGPU Advent Calender@21日目ですが、完全にネタ切れのため前回書き忘れていた部分を補足してお茶を濁します！
by jirohcl

まず、前回説明し忘れた重要なnoise関数。

float noise(float2 p,__global float2* rt,__global unsigned char* rti)
{
 float2 floor_p = floor(p);//float2(floor(p.x),floor(p.y))
 int i=(int)floor_p.x, j=(int)floor_p.y;
 float2 n00=rt[hash_idx(i,j,rti)];
 float2 n10=rt[hash_idx(i+1,j,rti)];
 float2 n01=rt[hash_idx(i,j+1,rti)];
 float2 n11=rt[hash_idx(i+1,j+1,rti)];

 float2 pf = p - floor_p;
 float sx= pow(pf.x,3)*(10-pf.x*(15-pf.x*6));
　　    float sy= pow(pf.y,3)*(10-pf.y*(15-pf.y*6));

 return bilerp(   dot(pf,n00)                   , dot(pf-(float2)(1.f,0.f),n10),
                         dot(pf-(float2)(0.f,1.f),n01) , dot(pf-(float2)(1.f,1.f),n11),
                         sx, sy);
}

となっており、パーティクルの座標と４点少しずらした座標でbilerpして求めます。
またランダムテーブル自体も

//host
float spin_t = 0.5f*0.03f/0.1f*time_;
for(unsigned int i=0; i<RANDOM_SIZE; ++i){
 float theta=rnd_table_spin_[i]*spin_t;
 float c=std::cos(theta),s=std::sin(theta);
 rnd_table_[i][0]= c*rnd_table_orig_[i][0]+s*rnd_table_orig_[i][1];
 rnd_table_[i][1]=-s*rnd_table_orig_[i][0]+c*rnd_table_orig_[i][1];
}

このように円を書くように更新していく事で、うずまくような動きになります。(つまりCurl)

float turb=g*d*0.0026f*noise((px-t*wind_velocity)/0.1f,rt,rti);
               ^^^^^^^

ノイズゲインを上げることによってこの動きは顕著になります。(前回の↑部分がゲイン)

といった所で前回の説明忘れ部分は終了です。
次にパフォーマンス比較を忘れていたので行いました。

100k particles
GPU          :   10ms
CPU single  : 183ms
CPU omp    :   93ms

CPU singleはシングルスレッド。CPU ompはOpenMPを用いた並列化の結果になります。
大体想定どおりの結果です。
CPU ompが若干遅いですが2coreのアレなCPUを使ってるので現行の4core辺りなら40-60ms位に落ち着くのではないでしょうか。

Curl Noise自体が並列化し易く、かつ入力データが多い例なのでOpenCLでの高速化が上手くはまるアルゴリズムでしたのでこのような(理想的な)結果に落ち着きましたが、何でもかんでもGPGPUで高速化できるものではありません。
例えば入力が多く、並列化が容易であっても分岐が多くなるとGPUでの高速化は上手くはまらないのでアルゴリズム審美眼を磨いて、「これはGPUでうまくいく！」といった勘を養うことがGPGPUでは重要ではないかなぁといった所で僕のノープランgdgd GPGPU Advent Calenderを終了しようと思います。
atnd立ち上げた山田さんに多謝！

Curl Noise with OpenCL

どうも再び@jirohclです。そしてGPGPU Advent Calender@12日目です！
人数が少なく、下手をすればもう一度まわってくるらしいのですが完全にネタも尽きてはたしてどうなる事やら。
まぁ愚痴はアンドロメダ星雲にでもすっ飛ばしておくとして、今回はこんなものを実装してみましたネタです。

Curl Noiseと聞いてピンと来る方も多いと思いますが、今をときめくスクウェアエニックス様が作ったAgni's Philosophyでも使われていたパーティクルの技法です。
主な特徴としては、パーティクル同士の判定を行わないために超高速、剛体との衝突がなんかそれっぽく見える、といった所でしょうか。

これをどうOpenCLに持っていくかというと、パーティクルの更新が都合よく並列化に向いていたためupdate関数をカーネル化してやりました。

__kernel void update(__global float2*        in,
		     __global float2*        rx,
		     __global float *         t,
		     __global float2*        rt,
		     __global unsigned char* rti){
	int const id = get_global_id(0);
	float2 mid = in[id];
	mid = mid+0.5f*curl_vel(in[id],rx,t[0],rt,rti)*0.04f;
	in[id] += curl_vel(mid,rx,t[0],rt,rti)*0.04f;
}

ここで特筆する点は乱数テーブル(rt)をホストから渡している所でしょうか。
この乱数テーブルは円状に配置された点をswapしてテーブル化してあります。

次に

float2 curl_vel(float2 x,__global float2* rx,float t,__global float2* rt, __global unsigned char* rti)
{
	float const delta_x = (1e-4);
	float2 v;
	v.x =  -(potential((float2)(x.x,x.y+delta_x),rx,t,rt,rti)
	       - potential((float2)(x.x,x.y-delta_x),rx,t,rt,rti))/(2*delta_x);
	v.y = (potential((float2)(x.x+delta_x,x.y),rx,t,rt,rti)
	     - potential((float2)(x.x-delta_x, x.y),rx,t,rt,rti))/(2*delta_x);
	return v;
}

速度の計算です、パーティクルの現在位置から微小距離ずらし、ポテンシャル関数を用いて速度を求めます。

float potential(float2 px,__global float2* rx,float t,__global float2* rt, __global unsigned char* rti)
{
	float2 wind_velocity = (float2)(0.1f,0.5f);
	float const inv_env_sqr = 1/sqrt(0.17f);
	float numer= inv_env_sqr*cross2(px,wind_velocity);
	float denom= inv_env_sqr;
        //---------------hard coded rigid X( ---
	float2 dx=px-rx[0];//centrer
	float psi=cross2(px,rx[1]);
	float phi=((sqrt(dx.x*dx.x)+sqrt(dx.y*dx.y))-0.06f);//radius;
	float m=1/(sqrt(phi)+1e-18);
	numer+=m*psi;
	denom+=m;
	//--------------------------------------
	float base_psi=numer/denom;

	float d=phi/0.1f;
	float g=smooth_step(1.1-px.x);
	float turb=g*d*0.0026f*noise((px-t*wind_velocity)/0.1f,rt,rti);
	return base_psi+turb;
}

そしてこのポテンシャル関数の中で、回転するノイズ関数、剛体や風等の周りからの力の付与を計算することで流体を表現するらしいです。

実行結果


実装に疲れたので説明はおざなりですが中々面白い技術なので皆様試してみてはいかがでしょうか？

Windows用の実行ファイルになり、かつ他の環境で動くかどうか怪しいですが一応デモも置いておきます。
※PCが壊れた等の責任は負いかねますので悪しからず

Hello Oldschool GPGPU World

お久しぶりです。@jirohclです。
特に書くことも思い浮かばず日が過ぎ、気がつけばまたAdvent Calenderの季節になりました。

今回、私が参加するのはこちら、GPGPU Advent Calender@3日目となります。
この記事は温故知新と聞こえの良い事を言いつつ、ネタが思い浮かばなかった苦肉の策、
タイトルの通り古臭いGPGPUをやってみようと思います。

古臭いGPGPUとは何かと言うと、プログラマブルシェーダ、特にピクセルシェーダを用いて汎用計算をやってみましょうというアレです。
これを知ることでCUDA/OpenCL等は裏側で一体何をしているのか押し知ることが出来、また効率的な実装の方法も見えてくるはず。（建前）

まず、GPGPUをするために何が必要かというとご存知の通り
1.デバイスへの入力
2.処理関数
3.ホストへの出力
これら３つ。

まず、1の入力はテクスチャを使うのがポピュラーな方法です。
次に２の処理関数に関してはフラグメントシェーダーで頑張ります。
３はglReadPixels等のレンダリング後の画像を読み込む関数を使うことで結果を取り出します。

ではまずデバイス側のコード(行列積)を見てみましょう。

//vs
uniform mat4 mvp;

in vec3 in_pos;
in vec2 in_tex;

out vec2 out_tex;

void main(void)
{
	gl_Position = mvp * vec4(gl_Vertex.xyz,0.f);
	out_tex = in_tex;
}

//fs
uniform int mat_size;

uniform sampler2D mat_1;
uniform sampler2D mat_2;

in vec2 texcoord;

void main (void)
{
	float result = 0.f;
	for(int i = 0;i<mat_size;++i)
	{
		result += texture2D( mat_1,vec2(i/(float)mat_size,texcoord.y) )
			* texture2D( mat_2,vec2(texcoord.y,i/(float)mat_size) );
	}
	gl_FragColor = vec4(result,1.f);
}

バーテクスシェーダ自体は何も特別なことは行わず、フラグメントシェーダで積を演算しています。
ここでフラグメントシェーダでの演算ですが実は一つ大きな問題があり、
フラグメントシェーダー内での処理には時間制限があるために行列のサイズが大きくなると計算が正常に終了しなくなります。OpenCLを使ったことがある方はきっと身に覚えがあるかと思います。アレです。

では次にホストコードを、と思いましたが全部羅列すると酷いことになるので一部

//データ作成
GLuint mat_lhs;
::glGenTextures( 1, &mat_lhs );
::glBindTexture( GL_TEXTURE_2D, mat_lhs );
::glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,mat_size,mat_size,0,GL_LUMINANCE,GL_FLOAT,data);
::glBindTexture( GL_TEXTURE_2D, 0 );

//実行
::glEnableClientState(GL_VERTEX_ARRAY);
::glEnableClientState(GL_TEXTURE_COORD_ARRAY);
::glBindBuffer(GL_ARRAY_BUFFER, texed_quad);
::glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, index_buffer);
::glDrawElements(GL_TRIANGLES,buffer_size/sizeof(GLuint),GL_UNSIGNED_INT,0);
::glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
::glBindBuffer(GL_ARRAY_BUFFER, 0);
::glDisableClientState(GL_VERTEX_ARRAY);
::glDisableClientState(GL_TEXTURE_COORD_ARRAY);
::glFlush();
::glReadPixels(0,0,mat_size,mat_size,GL_RGBA,GL_FLOAT,result);

デバイスへのデータ転送、カーネルの実行する”だけ”でもこれだけの決まりごとの命令が必要になります。
また、この他にもシェーダの初期化、ユニフォーム、イン/アウトへの関連付け、等々面倒くさい作業が目白押しとなっております。

CUDA/OpenCLが難しそうと思った方、古くさいGPGPUを動かすためのコードを見ると、現在の環境は大分めぐまれてると思いませんか？
思った方は是非GPGPU Advent Calenderに登録して記事を書いてみましょう！