[CUDA] Simulating Heat Transfer

이번 글에서는 CUDA를 이용하여 열전달 시뮬레이션을 병렬로 계산하는 예제에 대하여 다루도록 하겠다.

J. Sanders와 E. Kandrot의 CUDA by Example의 Chapter 7. Texture Memory 를 참고하여 재구성하였으며, 좀 더 심도있는 설명을 하기 위해 노력하였다.

[실행환경]

OS: Microsoft Windows 7 Professional 64

CUDA: 7.0

IDE: Microsoft Visual Studio 2012

Include: book.h, cpu_anim.h

book.h

cpu_anim.h

---

1.  Texture Memory Overview

Texture memory는 constant memory와 마찬가지로 온-칩(on-chip)에 캐시(cache)된다. 따라서 오프-칩(off-chip) DRAM으로의 메모리 요청을 줄여 고성능의 밴드 폭을 제공한다.

특히 texture cache는 메모리 액세스 패턴(memory access pattern)이 공간적으로 몰려있는 형태의 그래픽스 어플리케이션에 적합하도록 설계되었으며 아래 그림과 같이 각 쓰레드는 다른 쓰레드의 "근처"에 위치해 있다.

    [그림 1.] Texture memory에 대한 메모리 액세스 패턴

Texture memory의 CUDA 아키텍쳐 상의 구조적 위치는 다음 그림을 참고하도록 한다.

    [그림 2.] CUDA 아키텍쳐 상의 texture memory의 구조적 위치

[Texture Memory 읽기]

Texutre를 읽는 것을 texture fetching이라고 하며, 이에 대한 자세한 가이드는 CUDA Programming Guide B.8 Texture Functions를 참고한다.

[Texture Memory의 장점]

(1) 근방 데이터 값들의 선형 보간

(2) 자동 정규화

• 부호가 없는 경우: [0,+1]

• 부호가 있는 경우: [-1,+1]

(3) 자동 경계 처리

[Texture Memory의 단점]

Texture memory는 기본적으로 read-only이다. 그러나, 2D texture에 대해서는 메모리를 쓸 수 있는 방법이 있다.

[Texture Memory는 언제 사용되는가?]

(1) Memory access pattern이 공간 지역성(spatial locality)을 형성할 때

(2) 정확한 Read access pattern을 예측하기 어려울 때

(3) 그래픽스 파이프라인을 이용하여 데이터를 가시화할 때

[Texture Memory를 이용한 코딩 순서]

(1) CUDA에서 texture memory를 선언한다.

(2) CUDA에서 texture memory를 texture reference에 binding한다.

(3) CUDA kernel에서 texture reference로부터 texture memory를 읽는다.

(4) CUDA에서 texture reference부터 texture memory를 unbinding한다.

---

2.  Theoretical Background of Heat Transfer

열전달에 대한 간단한 이론적 배경을 살펴보고자 한다.

열전도에 대한 지배방정식(governing equation)은 Fourier의 열전도 법칙이며 다음과 같다:

여기서,

는 단위면적당 열유량 (heat flow, W/m²),

는 열전도율(thermal conductivity, W/m·k),

는 gradient operator,

는 온도(temperature, K)

이다.

부호 (-)는 열전도가 일어나는 방향성을 의미하는데, 열이 흐름이 온도가 높은데서 낮은데로 흐르기 때문이다.

2차원 문제에 대하여 이를 적용하면 다음과 같다:

Equation 1.

만약 heater에 의한 heat generation이 없는 steady-state에 대하여 continuity equation을 얻을 수 있다:

Equation 2.

여기서, 열전도율은 homogeneous하다고 가정하였다. 특히 Equation 2.는 potential problem의 일종으로 Laplace equation이라고 한다:

이 글에서는 steady-state의 Laplace equation을 푸는 목적은 아니며, transient problem에 대하여 물리적으로 엄밀하지는 않지만 CUDA에서의 texture memory를 활용하는 것에 초점을 맞추어 다음과 같은 수식으로 열전도를 시뮬레이션 하도록 하겠다:

Equation 3.

Equation 3.은 다음과 같이 풀어서 쓸 수 있다:

Equation 4.

---

3.  Source Code

(1) 상수 및 전역변수 선언

#define DIM 1024
#define PI 3.1415926535897932f
#define MAX_TEMP 1.0f
#define MIN_TEMP 0.0001f
#define SPEED   0.25f
#define ITERATIONS_PER_CALCULATION 90

// these exist on the GPU side
texture<float,2>  texConstSrc;
texture<float,2>  texIn;
texture<float,2>  texOut;

솔루션(픽셀) 공간의 해상도(DIM)는 1024*1024이며, SPEED는 열전도율로서 값이 클수록 열전도 속도가 빠르다.

texture<float, 2D>는 GPU에 float type으로 2차원 texture reference를 선언한다.

ITERATIONS_PER_CALCULATION 상수는 90번 반복을 계산 1 사이클로 정의한다는 의미이다.

(2) main 함수

DataBlock은 업데이트 루틴에 필요한 전역 구조체이다.

// globals needed by the update routine
struct DataBlock
{
    unsigned char   *output_bitmap;
    float           *dev_inSrc;
    float           *dev_outSrc;
    float           *dev_constSrc;
    CPUAnimBitmap   *bitmap;

    cudaEvent_t     start, stop;
    float           totalTime;
    float           frames;
};

DataBlock 구조체 인스턴스를 선언한다.

DataBlock   data;

이외에 다음과 같이 인스턴스들을 선언한다.

CPUAnimBitmap bitmap( DIM, DIM, &data );
data.bitmap = &bitmap;
data.totalTime = 0;
data.frames = 0;

bitmap을 1024×1024 크기로 생성한다.

GPU 계산의 성능을 측정하기 위해 CUDA Event에 다음과 같이 등록한다.

HANDLE_ERROR( cudaEventCreate( &data.start ) );
HANDLE_ERROR( cudaEventCreate( &data.stop ) );

변수 imageSize를 선언하고 GPU에 메모리를 생성한다.

int imageSize = bitmap.image_size();
HANDLE_ERROR( cudaMalloc( (void**)&data.output_bitmap, imageSize ) );

GPU에 다음 세 가지의 메모리를 생성한다.

// assume float == 4 chars in size (ie rgba)
HANDLE_ERROR( cudaMalloc( (void**)&data.dev_inSrc, imageSize ) );
HANDLE_ERROR( cudaMalloc( (void**)&data.dev_outSrc, imageSize ) );
HANDLE_ERROR( cudaMalloc( (void**)&data.dev_constSrc, imageSize ) );

Channel descriptor를 생성한다.

cudaChannelFormatDesc desc = cudaCreateChannelDesc<float>();

각 texture reference로 texture binding한다.

HANDLE_ERROR( cudaBindTexture2D( NULL, texConstSrc, data.dev_constSrc, desc, DIM, DIM,sizeof(float) * DIM ) );
HANDLE_ERROR( cudaBindTexture2D( NULL, texIn, data.dev_inSrc, desc, DIM, DIM, sizeof(float) * DIM ) );
HANDLE_ERROR( cudaBindTexture2D( NULL, texOut, data.dev_outSrc, desc, DIM, DIM, sizeof(float) * DIM ) );

cudaBindTexture2D의 프로토타입은 다음과 같다:

cudaError_t cudaBindTexture2D	(
	size_t * offset,
	const struct texture< T, dim, readMode > & tex,
	const void * devPtr,
	const struct cudaChannelFormatDesc & desc,
	size_t width,
	size_t height,
	size_t pitch
)

offset: byte 단위의 오프셋,

tex:  binding 할 texture reference,

devPrt:  device의 2D 메모리 영역,

desc:  channel 포맷,

width:  texture 단위의 넓이

height: texture 단위의 높이

pitch:  byte 단위의 pitch

Constant data를 초기화 한다. Constant data는 heater의 위치 범위와 그 범위에서의 heat generation 값을 지정한다. temp 포인터 변수에 heater의 범위와 heat generation 값을 정의한 후, cudaMemcpy 함수를 통해 device memory, 즉 GPU의 메모리 영역(data.dev_constSrc)으로 복사한다:

// initialize the constant data
float *temp = (float*)malloc( imageSize );

for( int i = 0 ; i < DIM*DIM ; i++ )
{
    temp[i] = 0;
    int x = i % DIM;
    int y = i / DIM;

	// sets heater
    if ( ( x > 300 ) && ( x < 600 ) && ( y > 310 ) && ( y < 601 ) )
        temp[i] = MAX_TEMP;
}

temp[DIM*100 + 100] = ( MAX_TEMP + MIN_TEMP ) / 2;
temp[DIM*700 + 100] = MIN_TEMP;
temp[DIM*300 + 300] = MIN_TEMP;
temp[DIM*200 + 700] = MIN_TEMP;

for( int y = 800 ; y < 900 ; y++ )
{
    for( int x = 400 ; x < 500 ; x++ )
	{
        temp[x + y*DIM] = MIN_TEMP;
    }
}

HANDLE_ERROR( cudaMemcpy( data.dev_constSrc, temp, imageSize, cudaMemcpyHostToDevice ) );

Input data를 초기화 한다. Input data는 heater 영역을 제외한 나머지 영역으로 모두 MIN_TEMP로 값이 정의된다. 마찬가지로 cudaMemcpy 함수를 통해 device memory, 즉 GPU의 메모리 영역(data.dev_inSrc)으로 복사한다:

// initialize the input data
for( int y = 800 ; y < DIM ; y++ )
{
    for( int x = 0 ; x < 200 ; x++ )
    {
        temp[x + y*DIM] = MAX_TEMP;
    }
}

HANDLE_ERROR( cudaMemcpy( data.dev_inSrc, temp, imageSize, cudaMemcpyHostToDevice ) );

임시 메모리 포인터 temp에 대한 메모리를 해제한다:

free( temp );

OpenGL을 이용한 animation loop을 정의한다:

bitmap.anim_and_exit( (void (*)(void*,int))anim_gpu, (void (*)(void*))anim_exit );

anim_and_exit()은 구조체 CPUAnimBitmap의 멤버 메써드이며, 헤더 cpu_anim.h에 정의되어 있다.

anim_and_exit()의 첫번째 인자는 함수 포인터와 정수형을 인자로 갖는 함수 포인터(anim_gpu())이며, 두번째 인자는 함수 포인터를 인자로 갖는 함수 포인터(anim_exit())이다.

cpu_anim.h의 anim_and_exit() 멤버 메써드를 살펴보면 다음과 같다:

void anim_and_exit( void (*f)(void*,int), void(*e)(void*) ) {
    CPUAnimBitmap**   bitmap = get_bitmap_ptr();
    *bitmap = this;
    fAnim = f;
    animExit = e;
    // a bug in the Windows GLUT implementation prevents us from
    // passing zero arguments to glutInit()
    int c=1;
    char* dummy = "";

    glutInit( &c, &dummy );
    glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGBA );
    glutInitWindowSize( width, height );
    glutCreateWindow( "bitmap" );
    glutKeyboardFunc(Key);
    glutDisplayFunc(Draw);

    if (clickDrag != NULL)
        glutMouseFunc( mouse_func );

    glutIdleFunc( idle_func );
    glutMainLoop();
}

fAnim과 animExit는 구조체 CPUAnimBitmap의 void 포인터 멤버이며, 함수 핸들러로 이해하면 된다. 구조체에 다음과 같이 선언되어 있다:

void (*fAnim)(void*,int);
void (*animExit)(void*);

glutIdleFunc() 함수는 인자로 idle_func를 넘겨주는데 이것은 일종의 callback 함수이며 다음과 같이 구현되어 있다:

// static method used for glut callbacks
static void idle_func( void ) {
    static int ticks = 1;
    CPUAnimBitmap*   bitmap = *(get_bitmap_ptr());
    bitmap->fAnim( bitmap->dataBlock, ticks++ );
    glutPostRedisplay();
}

anim_exit()는 Esc 키를 누르면 OpenGL main loop을 빠져나가기 위한 함수이다:

// static method used for glut callbacks
static void Key(unsigned char key, int x, int y) {
    switch (key) {
        case 27:
            CPUAnimBitmap*   bitmap = *(get_bitmap_ptr());
            bitmap->animExit( bitmap->dataBlock );

            //delete bitmap;
            exit(0);
    }
}

(3) anim_gpu 함수

/* @ function anim_gpu */
void anim_gpu( DataBlock *d, int ticks )
{
    HANDLE_ERROR( cudaEventRecord( d->start, 0 ) );
    dim3    blocks( DIM/16, DIM/16 );
    dim3    threads( 16, 16 );
    CPUAnimBitmap  *bitmap = d->bitmap;

    // since tex is global and bound, we have to use a flag to
    // select which is in/out per iteration
    volatile bool dstOut = true;

    for( int i = 0 ; i < ITERATIONS_PER_CALCULATION ; i++ )
	{
        float *in, *out;

        if( dstOut )
		{
            in  = d->dev_inSrc;
            out = d->dev_outSrc;
        }
		else
		{
            out = d->dev_inSrc;
            in  = d->dev_outSrc;
        }

        copy_const_kernel<<<blocks,threads>>>( in );
        blend_kernel<<<blocks,threads>>>( out, dstOut );
        dstOut = !dstOut;
    }

    float_to_color<<<blocks,threads>>>( d->output_bitmap, d->dev_inSrc );

    HANDLE_ERROR( cudaMemcpy( bitmap->get_ptr(),
                              d->output_bitmap,
                              bitmap->image_size(),
                              cudaMemcpyDeviceToHost ) );

    HANDLE_ERROR( cudaEventRecord( d->stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( d->stop ) );

    float   elapsedTime;

    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime, d->start, d->stop ) );

    d->totalTime += elapsedTime;
    ++d->frames;

    printf( "Average Time per frame:  %3.1f ms\n", d->totalTime/d->frames  );
}

Heat transfer simulation 코드는 매 프레임 마다 이 함수를 호출한다.

HANDLE_ERROR( cudaEventRecord( d->start, 0 ) );

에 의하여 함수 호출 시 성능 측정을 위해 먼저 Event Record를 시작한다. 그 다음 block과 thread의 크기를 설정하는데,

dim3    blocks( DIM/16, DIM/16 );

은 Grid 1개에 할당될 blocks 수를 지정한다. 지정된 blocks의 수는 DIM / 16 = 1024 / 16 = 16개이다.

dim3    threads( 16, 16 );

은 Block 1개에 할당될 threads 수를 지정한다.

dstOut은 일종의 flag인데 중요한 의미를 갖는다. 이 flag에 의하여 매번 계산할 때마다 in과 out이 교환된다. 이는 Equation 4.에서 해답을 찾을 수 있다.

즉, 는 와 로부터 계산되며, NEW는 새로 계산될 이미지(frame)을, OLD는 이전에 계산된 이미지(frame)이다. 이는 곧 최초에 out은 in으로부터 계산되며, 그 다음 계산에서는 in으로부터 out이 계산되는 방식이다.

그래서 for loop의 마지막 라인을 보면 dstOut의 상태를 계속하여 바꿔주는 것을 볼 수 있다:

dstOut = !dstOut;

다음은 CUDA kernel에서 constant 메모리를 복사하는 함수이다:

copy_const_kernel<<<blocks,threads>>>( in );

[copy_const_kernel 함수]

/* @ function copy_const_kernel */
__global__ void copy_const_kernel( float *iptr )
{
    // map from threadIdx/BlockIdx to pixel position
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    float c = tex2D( texConstSrc, x, y );

    if( c != 0 )
        iptr[offset] = c;
}

이 kernel 함수를 호출하는 함수는 anim_gpu() 함수이다:

copy_const_kernel<<<blocks,threads>>>( in );

copy_const_kernel() 함수는 texture reference texConstSrc로부터 지정된 위치(x, y)의 texture의 데이터를 얻어오고 이 데이터를 지정된 device memory에 저장한다. Texture 데이터를 얻어오는 함수의 프로토타입은 다음과 같다:

char tex2D(texture<char, cudaTextureType2D, cudaReadModeElementType> t, float x, float y)

[blend_kernel 함수]

anim_gpu() 함수는 GPU의 blend_kernel() 함수를 호출한다:

blend_kernel<<<blocks,threads>>>( out, dstOut );

blend_kernel() 함수는 다음과 같다:

/* @ function blend_kernel */
__global__ void blend_kernel( float *dst, bool dstOut )
{
    // map from threadIdx/BlockIdx to pixel position
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    float t, l, c, r, b;

    if( dstOut )
	{
        t = tex2D( texIn, x, y - 1 );
        l = tex2D( texIn, x - 1, y );
        c = tex2D( texIn, x, y );
        r = tex2D( texIn, x + 1, y );
        b = tex2D( texIn, x, y + 1 );
    } else
	{
        t = tex2D( texOut, x, y - 1 );
        l = tex2D( texOut, x - 1, y );
        c = tex2D( texOut, x, y );
        r = tex2D( texOut, x + 1, y );
        b = tex2D( texOut, x, y + 1 );
    }

    dst[offset] = c + SPEED * ( t + b + r + l - 4 * c );
}

Equation 4.를 계산하기 위하여 주어진 texel(texture element, 여기서는 하나의 picture element(pixel)로 간주한다)에 대하여 이웃의 texel 데이터를 얻어온다.

이 때 dstOut이 true일 경우에는 texIn texture reference을부터, dstOut이 false일 경우에는 texOut texture reference을부터 데이터를 얻어온다.

float으로 선언된 t, l, c, r, b는각각 주어진 texel을 기준으로 Top, Left, Center, Right, Bottom의 texel을 의미한다.

dst[offset] = c + SPEED * ( t + b + r + l - 4 * c );

은 Equation 4.를 계산하는 것이다.

다시 anim_gpu()로 돌아오자. for loop문을 통해 90번의 반복 계산이 완료되면 book.h 헤더 파일에 정의된 float_to_color() 함수로 이름 그대로 float을 color로 변환한다.

float_to_color<<<blocks,threads>>>( d->output_bitmap, d->dev_inSrc );

이 함수는 CUDA Toolkit에 정의된 함수가 아님을 유의하도록 한다. 함수의 프로토타입은 다음과 같다:

__global__ void float_to_color( unsigned char *optr, const float *outSrc )
{
    // map from threadIdx/BlockIdx to pixel position
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    float l = outSrc[offset];
    float s = 1;
    int h = (180 + (int)(360.0f * outSrc[offset])) % 360;
    float m1, m2;

    if (l <= 0.5f)
        m2 = l * (1 + s);
    else
        m2 = l + s - l * s;
    m1 = 2 * l - m2;

    optr[offset*4 + 0] = value( m1, m2, h+120 );
    optr[offset*4 + 1] = value( m1, m2, h );
    optr[offset*4 + 2] = value( m1, m2, h -120 );
    optr[offset*4 + 3] = 255;
}

코드에서 알 수 있듯이 float 포인터인 outSrc의 데이터로부터 R, G, B, A 채널값을 계산하여 optr에 저장한다.

이제 GPU device memory로부터 CPU host memory로 메모리를 복사한다:

HANDLE_ERROR( cudaMemcpy( bitmap->get_ptr(), d->output_bitmap, bitmap->image_size(), cudaMemcpyDeviceToHost ) );

그 다음 현재 image frame에 대한 CUDA Event를 종료한다. (1개의 image frame 계산을 위해 90번의 iteration이 수행되었다.)

HANDLE_ERROR( cudaEventRecord( d->stop, 0 ) );

그리고나서 가장 최근의 cudaEventRecord() 호출에 대한 모든 device 작업이 마칠 때까지 대기한다.

HANDLE_ERROR( cudaEventSynchronize( d->stop ) );

마지막으로 해당 image frame 계산이 걸린 시간을 체크하고, frame count를 업데이트 한다:

float   elapsedTime;

HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime, d->start, d->stop ) );

d->totalTime += elapsedTime;

++d->frames;

printf( "Average Time per frame:  %3.1f ms\n", d->totalTime/d->frames  );

---

4.  Result