目次 >> CUDA >> 拡散方程式

CUDA - 拡散方程式を解く

ここでは、CUDAを使って拡散方程式を解いてみる。
解き方はもっとも単純な差分法を使っている。

拡散方程式を解く（1次元）

まずは単純に、diff1d.cuというテキストファイルに、

#include <stdio.h>
#include <cutil.h>

#include "diff1d_kernel.cu"

int main( int argc, char** argv) 
{
    //デバイスの初期化
    CUT_DEVICE_INIT(argc, argv);

    //結果書き込み用ファイルのオープン
    FILE *fp=fopen("result.txt","w");
    //タイマーを作成して計測開始
    unsigned int timer = 0;
    CUT_SAFE_CALL( cutCreateTimer( &timer));
    CUT_SAFE_CALL( cutStartTimer( timer));

    //メインメモリ上にfloat型のデータを100個生成する
    float* h_idata = (float*) malloc(sizeof( float) * 100);
    //初期条件をセット
    for( int i = 0; i < 50; i++) 
        h_idata[i] = 0;
    for( int i = 50; i < 100; i++) 
        h_idata[i] = 1;

        //デバイス上（ビデオカードのこと）にも同じくfloat型100個分のメモリを確保する
    float* d_idata;
    CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata, sizeof( float) * 100 ));

    dim3  grid( 1, 1, 1);
    //100は100個並列であるため
    dim3  threads(100, 1, 1);
        
    //メインメモリからデバイスのメモリにデータを転送する
    CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, sizeof( float) * 100 , cudaMemcpyHostToDevice) );

    for (int t=0;t<100;t++)
    {
        //ここでGPUを使った計算が行われる
        diff1dKernel<<< grid, threads, sizeof( float) * 200 >>>( d_idata);
                
        //デバイスからメインメモリ上に実行結果をコピー
        CUDA_SAFE_CALL( cudaMemcpy( h_idata, d_idata, sizeof( float) * 100, cudaMemcpyDeviceToHost) );

        //実行結果を表示
        for (int i=0;i<100;i++)
        {
            fprintf(fp,"%f\t",h_idata[i]);
        }
        fprintf(fp,"\n");
    }
    //タイマーを停止しかかった時間を表示
    CUT_SAFE_CALL( cutStopTimer( timer));
    printf("Processing time: %f (ms)\n", cutGetTimerValue( timer));
    CUT_SAFE_CALL( cutDeleteTimer( timer));
    //各種メモリを解放
    free( h_idata);
    CUDA_SAFE_CALL(cudaFree(d_idata));
    fclose(fp);
    //終了処理
    CUT_EXIT(argc, argv);
    return 0;
}

diff1d_kernel.cuの方はというと、

__global__ void diff1dKernel( float* g_idata) 
{
    extern __shared__ char sharedmem[]; 
    float * sdata = (float *) sharedmem; 
    float * sdata2 = (float *) (sharedmem+sizeof(float)*100);

    // スレッドIDを取得
    const unsigned int tid = threadIdx.x;

    //グローバルメモリから入力データの読み込み
    sdata[tid] = g_idata[tid];
    __syncthreads();

    //ここで計算を行う
    const float Dfu=1;
    const float dt=0.2;
    const float dx=1;
    float Dfudtdx2=Dfu*dt/(dx*dx);
    for (int n=0;n<100;n++)
    {
        if (tid==0)
            sdata2[tid]=sdata[tid]+(2*sdata[tid+1]-2*sdata[tid])*Dfudtdx2;
        else if (tid==99)
            sdata2[tid]=sdata[tid]+(2*sdata[tid-1]-2*sdata[tid])*Dfudtdx2;
        else
            sdata2[tid]=sdata[tid]+(sdata[tid-1]+sdata[tid+1]-2*sdata[tid])*Dfudtdx2;
        __syncthreads();
        sdata[tid]=sdata2[tid];
    }

    //グローバルメモリに結果を書き込む
    g_idata[tid] = sdata[tid];
}

これを、

nvcc diff1d.cu -lcutil32 -lkernel32 -o diff1d.exe

とコンパイルし、できたdiff1d.exeを実行する。

>diff1d.exe
Processing time: 45.219727 (ms)

Press ENTER to exit...

Mac OS Xでコンパイルする場合は、

$ nvcc -L/Developer/CUDA/lib diff1d.cu -lcutil -o diff1d

のようにコンパイルする。

拡散方程式を解く（2次元）

拡散方程式を2次元で解く場合、1次元の時とは違い系が大きいので1次元の時ような形でのShared Memoryは使えない（Shared Memoryは各ブロックあたり16kBしかないので、floatの場合で、4096個まで使用できる）。
また、前回は、ブロック数は1であったので、GPUのプロセッサをすべて使い切っていない。ここでは、複数のブロックを使ってより並列度を上げてみる。
ブロックを複数使った場合の問題は、カーネル内でブロック間の同期をとる方法が存在しない点である。そのため、下記のプログラムでは、1回計算するたびに、カーネルを終了し同期をとっている。

まず、diff2d.cuであるが、

#include <stdio.h>
#include <cutil.h>

#include "diff2d_kernel.cu"

//系のサイズ 100 X 100のグリッド上で拡散方程式を解く
const int X=100;
const int Y=100;

int main( int argc, char** argv) 
{
    //デバイスの初期化
    CUT_DEVICE_INIT(argc, argv);

    //結果書き込み用ファイルのオープン
    FILE *fp=fopen("result.txt","w");

    //タイマーを作成して計測開始
    unsigned int timer = 0;
    CUT_SAFE_CALL( cutCreateTimer( &timer));
    CUT_SAFE_CALL( cutStartTimer( timer));

    //メインメモリ上にfloat型のデータをX*Y個生成する
    float* h_idata = (float*) malloc(sizeof( float) * X*Y);
    //初期条件をセット
    for( int i = 0; i < X; i++) 
        for( int j = 0; j < X; j++)
            if((i-X/2)*(i-X/2)+(j-Y/2)*(j-Y/2)<10*10)
                h_idata[i*Y+j] = 1;
            else
                h_idata[i*Y+j] = 0;

    //デバイス上（ビデオカードのこと）にも同じくfloat型X*Y個分のメモリを確保する
    float* d_idata;
    CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata, sizeof( float) * X*Y ));
    //デバイス上（ビデオカードのこと）にfloat型X*Y個分の作業用メモリを確保する
    float* d_idata2;
    CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata2, sizeof( float) * X*Y ));

    //ブロック数を増やして並列度を上げる
    dim3  grid( 16, 1, 1);
    dim3  threads(256, 1, 1);
    
    //メインメモリからデバイスのメモリにデータを転送する
    CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, sizeof( float) * X*Y , cudaMemcpyHostToDevice) );

    for (int t=0;t<100;t++)
    {
        for (int n=0;n<10;n++)
        {
            //ここでGPUを使った計算が行われる
            diff2dKernel<<< grid, threads>>>( d_idata, d_idata2,X,Y);
            //作業用領域から書き戻す
            CUDA_SAFE_CALL( cudaMemcpy( d_idata, d_idata2, sizeof( float) * X*Y, cudaMemcpyDeviceToDevice) );
        }
        //デバイスからメインメモリ上に実行結果をコピー
        CUDA_SAFE_CALL( cudaMemcpy( h_idata, d_idata, sizeof( float) * X*Y, cudaMemcpyDeviceToHost) );
        //実行結果を表示
        for (int i=0;i<X;i+=2)
        {
            for (int j=0;j<Y;j+=2)
            {
                fprintf(fp,"%f\t",h_idata[i*Y+j]);
            }
        }
        fprintf(fp,"\n");
    }
    //タイマーを停止しかかった時間を表示
    CUT_SAFE_CALL( cutStopTimer( timer));
    printf("Processing time: %f (ms)\n", cutGetTimerValue( timer));
    CUT_SAFE_CALL( cutDeleteTimer( timer));
    //各種メモリを解放
    free( h_idata);
    free( h_idata2);
    CUDA_SAFE_CALL(cudaFree(d_idata));
    fclose(fp);
    //終了処理
    CUT_EXIT(argc, argv);
    return 0;
}

diff2d_kernel.cuは、

__global__ void diff2dKernel(float* g_idata, float* g_idata2,int X, int Y)
{
    const unsigned int tid = threadIdx.x;
    const unsigned int bid = blockIdx.x;
    const unsigned int bdim = blockDim.x;
    const unsigned int gdim = gridDim.x;
    int step=bdim*gdim;
    int num=X*Y;

    //ここで計算を行う
    const float Dfu=1;
    const float dt=0.2;
    const float dx=1;
    float Dfudtdx2=Dfu*dt/(dx*dx);
    for (int id=bid * bdim + tid;id<num;id+=step)
    {
        //境界は、nonflux境界条件
        if ((id)==0)//四隅
            g_idata2[id]=g_idata[id]+(g_idata[id+1]+g_idata[id+1]+g_idata[id+Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2;
        else if ((id)==Y-1)//四隅
            g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id-1]+g_idata[id+Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2;
        else if ((id)==X*Y-Y)//四隅
            g_idata2[id]=g_idata[id]+(g_idata[id+1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id-Y]-4*g_idata[id])*Dfudtdx2;
        else if ((id)==X*Y-1)//四隅
            g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id-1]+g_idata[id-Y]+g_idata[id-Y]-4*g_idata[id])*Dfudtdx2;
        else if ((id)<Y)//四辺
            g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id+1]+g_idata[id+Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2;
        else if((id)>X*Y-Y)//四辺
            g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id-Y]-4*g_idata[id])*Dfudtdx2;
        else if ((id)%Y==0)//四辺
            g_idata2[id]=g_idata[id]+(g_idata[id+1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2;
        else if ((id)%Y==Y-1)//四辺
            g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id-1]+g_idata[id-Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2;
        else
            g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2;
    }
}

上記のプログラムを下記のようにしてコンパイルする。

>nvcc diff2d.cu -lcutil32 -lcuda -o diff2d.exe

そして実行する。当方の環境では、下記のような実行時間となった。

>diff2d.exe
Processing time: 762.077759 (ms)

Press ENTER to exit...

結果は下記の通り。

この結果の表示にはMatlabを使った。一応そのコード。

data=load('result.txt');
for lp=1:9
subplot(3,3,lp)
z=data(lp*5,:);
z=reshape(z,50,50);
surf(z)
shading interp
daspect([50 50 1])
axis([1,50,1,50,0,1])
caxis([0 1])
end

本文中のFC4はFedora ProjectのFedora Core 4を、FC5はFedora Core 5を、FC6はFedora Core 6をopenSUSEはNovellのSUSE Linux OSSを表します。Fedora7以降は、単にFedora7、Fedora8、Fedora9、Fedora10、Fedora11、Fedora12、Fedora13、Fedora14、Fedora15と表示しています。Ubuntuは、必要に応じて10.04、11.04のようにバージョン番号をつけて区別しています。MandrivaはMandriva Linuxを表します。