GPUの演算能力を、物理シミュレーションとグラフィックスの融合に活用してみました。 本プロジェクトでは、流体解析の手法である「格子ボルツマン法（LBM: Lattice Boltzmann Method）」をCUDAで実装し、OpenGLを用いてその動態をリアルタイムに可視化するシミュレーターを作成しました。
中心には複雑な干渉を生む「星型」の障害物を配置。マウスドラッグによる滑らかな波紋の生成や、障害物による波の回折・干渉を、512×512のグリッド上で高フレームレートにシミュレートします。物理演算と描画の両方をGPU内で完結させることで、CPUでは到底不可能なリアルタイム性を実現しています。
技術詳細解説
本プログラムは、大きく分けて**「物理演算（LBM）」「パルス注入（入力）」「可視化」**の3つのフェーズで構成されています。
1. 格子ボルツマン法 (D2Q9モデル)
流体を微視的な粒子の集まりとして捉え、その速度分布関数の変化を計算します。
-
D2Q9モデル: 各格子点は9つの速度方向（静止、上下左右、斜め）を持ちます。
-
ストリーミング（移動）: 隣接する格子から粒子が移動してくるステップです。本コードでは、境界条件（星型の壁）に衝突した粒子を跳ね返す「Bounce-back法」を実装しています。
-
衝突（緩和）: 粒子同士の衝突により、局所的な平衡状態へ近づくステップです。定数
omega（緩和時間）が流体の粘性を決定します。
2. 滑らかなマウス入力（線形補間パルス）
マウスドラッグの軌跡を滑らかに表現するため、単なる点ではなく「線分」としてエネルギーを注入しています。
-
線分最短距離計算: 前フレームと現フレームの座標を結ぶ線分から、各格子点への最短距離を計算します。
-
ガウス分布による注入: 距離に応じた重み付け（
expf）を行い、密度を加算することで、航跡のような自然な波紋を生成します。
3. GPUダイレクト可視化
演算結果を一度もCPUに戻すことなく、GPU内で直接映像へ変換しています。
-
物理量のカラーマッピング: 格子点ごとの流速（
ux,uy）と密度（rho）の変化を色相に変換。tanhf関数を用いることで、激しい変化もサチュレーションを起こさず美しく描写します。 -
テクスチャ共有: CUDAで計算した画像バッファをOpenGLのテクスチャとしてバインドし、1枚の板ポリゴンに貼り付けて表示しています。
12コンパイルコマンドnvcc lbm_cuda_opengl.cu -o lbm_sim -arch=sm_61 -lGLEW -lGLU -lGL -lglut -lm123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154#include <GL/glew.h>#include <GL/freeglut.h>#include <cuda_runtime.h>#include <vector>#include <math.h>#include <stdio.h>const int WIDTH = 512, HEIGHT = 512, Q = 9;int win_w = WIDTH, win_h = HEIGHT;const float omega = 0.7f;bool is_dragging = false;int prev_mx = -1, prev_my = -1, cur_mx = -1, cur_my = -1;float *d_f1, *d_f2, *d_ux, *d_uy, *d_rho;int *d_boundary;uchar4 *d_img;std::vector<uchar4> h_output_buffer(WIDTH * HEIGHT);GLuint tex;__device__ const int d_cx[9] = {0, 1, 0, -1, 0, 1, -1, -1, 1};__device__ const int d_cy[9] = {0, 0, 1, 0, -1, 1, 1, -1, -1};__device__ const float d_w[9] = {4/9.f, 1/9.f, 1/9.f, 1/9.f, 1/9.f, 1/36.f, 1/36.f, 1/36.f, 1/36.f};const float h_w[9] = {4/9.f, 1/9.f, 1/9.f, 1/9.f, 1/9.f, 1/36.f, 1/36.f, 1/36.f, 1/36.f};// --- パルス生成：線形補間 ---__global__ void pulse_line_kernel(float *f, int x1, int y1, int x2, int y2, float mag) {int x = blockIdx.x * blockDim.x + threadIdx.x, y = blockIdx.y * blockDim.y + threadIdx.y;if (x >= WIDTH || y >= HEIGHT) return;float dx = (float)(x2 - x1), dy = (float)(y2 - y1);float t = fmaxf(0.0f, fminf(1.0f, ((x - x1) * dx + (y - y1) * dy) / (dx * dx + dy * dy + 1e-8f)));float px = x1 + t * dx, py = y1 + t * dy;float dist2 = (x - px) * (x - px) + (y - py) * (y - py);float weight = expf(-dist2 / 200.0f);if (weight > 0.01f) {for (int q = 0; q < Q; q++) f[q * WIDTH * HEIGHT + y * WIDTH + x] += d_w[q] * weight * mag;}}// --- LBM演算 ---__global__ void lbm_kernel(float *fin, float *fout, float *rho, float *ux, float *uy, const int *bound) {int x = blockIdx.x * blockDim.x + threadIdx.x, y = blockIdx.y * blockDim.y + threadIdx.y;if (x >= WIDTH || y >= HEIGHT) return;int idx = y * WIDTH + x;if (bound[idx]) return;float f[9], r = 0, vx = 0, vy = 0;int opp[9] = {0, 3, 4, 1, 2, 7, 8, 5, 6};for (int q = 0; q < Q; q++) {int sx = (x - d_cx[q] + WIDTH) % WIDTH, sy = (y - d_cy[q] + HEIGHT) % HEIGHT;f[q] = bound[sy * WIDTH + sx] ? fin[opp[q] * WIDTH * HEIGHT + idx] : fin[q * WIDTH * HEIGHT + sy * WIDTH + sx];r += f[q]; vx += d_cx[q] * f[q]; vy += d_cy[q] * f[q];}rho[idx] = r; ux[idx] = vx / (r + 1e-8f); uy[idx] = vy / (r + 1e-8f);for (int q = 0; q < Q; q++) {float cu = 3.f * (d_cx[q] * ux[idx] + d_cy[q] * uy[idx]);float feq = d_w[q] * r * (1.f + cu + 0.5f * cu * cu - 1.5f * (ux[idx] * ux[idx] + uy[idx] * uy[idx]));fout[q * WIDTH * HEIGHT + idx] = f[q] + omega * (feq - f[q]);}}// --- 可視化（障害物を灰色に設定） ---__global__ void visual_kernel(uchar4 *out, const float *ux, const float *uy, const float *rho, const int *bound) {int x = blockIdx.x * blockDim.x + threadIdx.x, y = blockIdx.y * blockDim.y + threadIdx.y, idx = y * WIDTH + x;if (x >= WIDTH || y >= HEIGHT) return;if (bound[idx]) {out[idx] = make_uchar4(100, 100, 100, 255); // 灰色に変更return;}float i = tanhf(sqrtf(ux[idx]*ux[idx] + uy[idx]*uy[idx]) * 30.f + fabsf(rho[idx] - 1.f) * 6.f);out[idx] = make_uchar4(i * 255, i * 130, 210 * (1 - i) + 40, 255);}void reset() {std::vector<float> h_f(Q * WIDTH * HEIGHT);for (int i = 0; i < WIDTH * HEIGHT; i++) for (int q = 0; q < Q; q++) h_f[q * WIDTH * HEIGHT + i] = h_w[q];cudaMemcpy(d_f1, h_f.data(), Q * WIDTH * HEIGHT * sizeof(float), cudaMemcpyHostToDevice);}void disp() {dim3 b(16, 16), g((WIDTH + 15) / 16, (HEIGHT + 15) / 16);if (is_dragging && cur_mx != -1) {pulse_line_kernel<<<g, b>>>(d_f1, prev_mx, prev_my, cur_mx, cur_my, 0.08f);prev_mx = cur_mx; prev_my = cur_my;}lbm_kernel<<<g, b>>>(d_f1, d_f2, d_rho, d_ux, d_uy, d_boundary);std::swap(d_f1, d_f2);visual_kernel<<<g, b>>>(d_img, d_ux, d_uy, d_rho, d_boundary);cudaMemcpy(h_output_buffer.data(), d_img, WIDTH * HEIGHT * sizeof(uchar4), cudaMemcpyDeviceToHost);glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, WIDTH, HEIGHT, GL_RGBA, GL_UNSIGNED_BYTE, h_output_buffer.data());glClear(GL_COLOR_BUFFER_BIT);glBegin(GL_QUADS); glTexCoord2f(0,0); glVertex2f(-1,-1); glTexCoord2f(1,0); glVertex2f(1,-1);glTexCoord2f(1,1); glVertex2f(1,1); glTexCoord2f(0,1); glVertex2f(-1,1); glEnd();glutSwapBuffers(); glutPostRedisplay();}void mouse(int b, int s, int x, int y) {if (b == GLUT_LEFT_BUTTON) {if (s == GLUT_DOWN) {is_dragging = true;cur_mx = prev_mx = x * WIDTH / win_w;cur_my = prev_my = (win_h - y) * HEIGHT / win_h;pulse_line_kernel<<<dim3((WIDTH+15)/16,(HEIGHT+15)/16), dim3(16,16)>>>(d_f1, cur_mx, cur_my, cur_mx, cur_my, 0.2f);} else is_dragging = false;}}void motion(int x, int y) {if (is_dragging) { cur_mx = x * WIDTH / win_w; cur_my = (win_h - y) * HEIGHT / win_h; }}// 星型判定関数bool is_star(float x, float y, float cx, float cy, float r, float inner_r) {float dx = x - cx, dy = y - cy;float dist = sqrtf(dx*dx + dy*dy);float angle = atan2f(dy, dx);// 5つの頂点を持つ星の形状パラメータfloat s = (angle + M_PI/2.0f) * 5.0f / (2.0f * M_PI);float f = s - floorf(s);float p = (f < 0.5f) ? f : 1.0f - f;float star_r = r * (1.0f - p * (r - inner_r) / r * 4.0f);return dist < star_r;}int main(int argc, char **argv) {glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);glutInitWindowSize(WIDTH, HEIGHT); glutCreateWindow("LBM: Star Obstacle");glewInit();cudaMalloc(&d_f1, Q*WIDTH*HEIGHT*4); cudaMalloc(&d_f2, Q*WIDTH*HEIGHT*4);cudaMalloc(&d_ux, WIDTH*HEIGHT*4); cudaMalloc(&d_uy, WIDTH*HEIGHT*4);cudaMalloc(&d_rho, WIDTH*HEIGHT*4); cudaMalloc(&d_boundary, WIDTH*HEIGHT*4);cudaMalloc(&d_img, WIDTH*HEIGHT*4);std::vector<int> h_b(WIDTH * HEIGHT, 0);for (int i = 0; i < WIDTH * HEIGHT; i++) {int x = i % WIDTH, y = i / WIDTH;// 外壁if (x == 0 || x == WIDTH - 1 || y == 0 || y == HEIGHT - 1) h_b[i] = 1;// 星型障害物 (中心 256, 256 付近、半径 60, 内径 25)if (is_star((float)x, (float)y, 256.0f, 256.0f, 65.0f, 25.0f)) h_b[i] = 1;}cudaMemcpy(d_boundary, h_b.data(), WIDTH * HEIGHT * 4, cudaMemcpyHostToDevice);reset();glGenTextures(1, &tex); glBindTexture(GL_TEXTURE_2D, tex);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, WIDTH, HEIGHT, 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);glEnable(GL_TEXTURE_2D);glutDisplayFunc(disp); glutReshapeFunc([](int w, int h) { glViewport(0,0,w,h); win_w=w; win_h=h; });glutMouseFunc(mouse); glutMotionFunc(motion);glutKeyboardFunc([](unsigned char k, int x, int y) { if (k == 'r' || k == ' ' || k == 'R') reset(); });glutMainLoop(); return 0;}This project is a real-time fluid dynamics simulator that implements the Lattice Boltzmann Method (LBM) by leveraging the parallel processing power of the GPU. By tightly integrating CUDA for physical computations and OpenGL for rendering, it achieves exceptionally smooth visualization of fluid behavior on a high-resolution grid.
Technical Highlights
-
Star-Shaped Obstacle: A gray, star-shaped obstacle is positioned at the center of the domain. Utilizing a mathematical Signed Distance Function (SDF) approach based on polar coordinates, the simulator accurately handles wave diffraction and interference against complex geometries.
-
Linearly Interpolated Drag Input: To ensure seamless wave generation regardless of mouse movement speed, a custom kernel interpolates between previous and current mouse coordinates, injecting energy along the continuous path.
-
Zero-Copy Visualization: Physical quantities are converted directly into pixel colors within GPU memory. These are then rendered as OpenGL textures without being transferred back to the CPU, eliminating data bottlenecks and ensuring high frame rates.
How It Works
-
Collision & Streaming: The particle distribution functions at each grid point move to neighboring nodes and collide to approach a local equilibrium state.
-
Boundary Conditions: Particles hitting the star-shaped wall undergo a “bounce-back” process, creating intricate wave patterns around the object.
-
Visualization: Changes in local velocity and density are mapped to a vibrant color palette using a
tanhactivation function to prevent color saturation.
-