简单的传热模型的动画模拟

/*
简单的二维热传导模拟：
首先，假设有一个矩形房间，将其分成一个网格，在网格中随机散布一些热源，
他们有着不同的固定温度，在给定了矩形网格及热源分布后可以计算网格中每
个单元的温度随时间的变化情况，为了简单热源单元本身的温度将保持不变，
在时间递进的每个步骤中，假设热量在某个单干及其邻接单元之间流动，如果
某个单元的临接单元的温度更高，那么热量将从临接单元传导到该单元，相反的
如果某个单元比临接单元温度高，那么它将变冷。
温度更新的计算步骤：
1）给定一个包含初始输入温度的网格，将其中作为热源的单元温度值复制到网格中
相应的单元中，这将覆盖这些单元之前计算的温度，因些也就确保了加热单元将保持
恒温这个条件。这个复制操作是在copy_const_kernel（）中执行的。
2）给定一个温度输入网格，根据公式7.2中的更新公式计算输出温度风格，这个更新
操作是在blend_kernel()中执行的。
3）将输入温度网格和输出温度网格交换，为下一个步骤的计算作好准备，当模拟下一
个时间步骤时，在步骤2）中计算得到的输出温度网格将成为步骤1）中的输入温度网格

*/
#include "cuda.h"
#include "book.h"
#include "cpu_anim.h"

#define DIM 1024
#define PI 3.1415926535f
#define MAX_TEMP 1.0f   //最高温度
#define MIN_TEMP 0.0001f    //最低温度
#define SPEED 0.25f  //热量的流动速度

//将作为热源的单元温度值复制到网格中相应的单元中，
__global__ void copy_const_kernel(float *iptr,const float *cptr){
 //
 int x=threadIdx.x+blockIdx.x*blockDim.x;
 int y=threadIdx.y+blockIdx.y*blockDim.y;
 int offset=x+y*blockDim.x*gridDim.x;

 if(cptr[offset]!=0)
  iptr[offset]=cptr[offset]; //把cptr的热源温度复制到iptr的输入单元中
}
//每个线程都负责计算一个单元，每个线程都将读取对应单元及其邻接单元的温度值
//执行前面给出的更新运算，即用上，下，左，右相邻的单元温度计算当前单元的温度
__global__ void blend_kernel(float *outSrc,const float *inSrc){
 //
 int x=threadIdx.x+blockIdx.x*blockDim.x;
 int y=threadIdx.y+blockIdx.y*blockDim.y;
 int offset=x+y*blockDim.x*gridDim.x;

 int left=offset-1;
 int right=offset+1;
 if(x==0) left++;  //左边界处理
 if(x==DIM-1) right--;   //右边界处理

 int top=offset-DIM;
 int bottom=offset+DIM;
 if(y==0) top+=DIM;    //上边界处理
 if(y==DIM-1) bottom-=DIM;   //下边界处理
//用四邻接单元计算当前单元的温度
 outSrc[offset]=inSrc[offset]+SPEED*(inSrc[top]+
             inSrc[bottom]+inSrc[left]+inSrc[right]-
       inSrc[offset]*4);
}

//数据块的定义
struct DataBlock{
 unsigned char *output_bitmap;  //输出位图数据的定义
 float *dev_inSrc;
 float *dev_outSrc;
 float *dev_constSrc;
 CPUAnimBitmap *bitmap;  //输出位图的定义

 cudaEvent_t start,stop;   //定义记录GPU处理时间的事件
 float totalTime;
 float frames;
};
//
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;

 for(int i=0;i<90;i++)
 {   //计算模拟过程的一个时间步
  copy_const_kernel<<<blocks,threads>>>(d->dev_inSrc,d->dev_constSrc);  //算法第一步：复制温度
  blend_kernel<<<blocks,threads>>>(d->dev_outSrc,d->dev_inSrc);   //算法第二步：计算输出温度
  swap(d->dev_inSrc,d->dev_outSrc);   //算法第三步，新的输出作为下次循环的输入
 }
 //float_to_color在本书带的库book.h中有相应的定义：功能是将温度与颜色建立映射关系
 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);
}
//释放显存空间
void anim_exit(DataBlock *d)
{
 cudaFree(d->dev_inSrc);
 cudaFree(d->dev_outSrc);
 cudaFree(d->dev_constSrc);

 HANDLE_ERROR(cudaEventDestroy(d->start));
 HANDLE_ERROR(cudaEventDestroy(d->stop));
}

int main(void)
{
 DataBlock data;
 CPUAnimBitmap bitmap(DIM,DIM,&data);
 data.bitmap=&bitmap;
 data.totalTime=0;
 data.frames=0;

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

 HANDLE_ERROR(cudaMalloc((void**)&data.output_bitmap,bitmap.image_size()));
 //
 HANDLE_ERROR(cudaMalloc((void**)&data.dev_inSrc,bitmap.image_size()));
 HANDLE_ERROR(cudaMalloc((void**)&data.dev_outSrc,bitmap.image_size()));
 HANDLE_ERROR(cudaMalloc((void**)&data.dev_constSrc,bitmap.image_size()));

 float *temp=(float*)malloc(bitmap.image_size());
 for(int i=0;i<DIM*DIM;i++)
 {
  temp[i]=0;
  int x=i%DIM;
  int y=i/DIM;
  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+100]=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,bitmap.image_size(),cudaMemcpyHostToDevice));
 
 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,bitmap.image_size(),cudaMemcpyHostToDevice));

 free(temp);
 //传递函数指针调用核函数处理动画数据
 bitmap.anim_and_exit((void(*)(void*,int))anim_gpu,(void(*)(void*))anim_exit);
}