基于 格子Boltzman方法的poiseuille流模拟

  

     跟着老师在做LBM流体模拟，感觉理论的东西特别的难，看了大半天理论，虽然还懂为什么要这么做，但是还是先用代码实现起来先，对比CPU和GPU的实现做一个比较。

    现在做的的第一个例子是poiseuille流，先用CPU做，才用GPU来做。

    先上题目：

    两平行板间的二维流场只存在沿x方向的流动，假设流动为定常且忽略体积力，见下图。

     

    格子Boltzman 方法（35*101，上下固壁上采用半程反弹，两边采样周期性边界）的图片：

     

     初态每个点的水平和垂直方向的速度都是0，每次碰撞后在水平方向上加了一个体力驱动，所以在水平方向上的的速度矢量图：

     

    中间的速度是最大的，上下两边的速度最小，形成了一个抛物线的形状，最后基本趋于一个稳态。

   

    CPU实验环境：VS2008、OpenGL库（已经在项目文件中）

    下面是代码实现

#include <stdio.h>#include <time.h>#include <iostream>#include "../common/cpu_anim.h"using namespace std;#define Width 1001#define Height 1001#define ImageScale 2          //图片宽度扩宽倍数#define ImageScaleU_x 20000   //水平速度放大倍数#define remainderX(x) ((x + Width) % Width)    //周期性边界#define remainderY(x) ((x + Height) % Height)  //周期性边界#define W0 0.4444444444444444444444444 //(4.0/9.0)#define W1 0.1111111111111111111111111 //(1.0/9.0)#define W2 0.02777777777777777777777778 //(1.0/36.0)double ForceX=1E-6;double ForceY=0;double Gamma[9] ={ 0, 1.0/3, 1.0/3, 1.0/3, 1.0/3, 1.0/12, 1.0/12, 1.0/12, 1.0/12 };double dAlpha[]={W0,W1,W1,W1,W1,W2,W2,W2,W2};int ix[9]={0, 0,1,0,-1, 1,1,-1,-1};int iy[9]={0,-1,0,1, 0,-1,1, 1,-1};double Den=1.0;double dRevTau=1.0/0.75;       //弛豫时间int inverse[9]={0,3,4,1,2,7,8,5,6};//相反方向struct DataBlock{float *m_df;  //设备上每个点的数据float *u_x;     //host上x方向上的宏观速度float *u_y;     //host上y方向上的宏观速度float *dDen;    //每个点的密度CPUAnimBitmap *bitmap;};double feq(int Order,double Den,double Velx,double Vely)  //泊肃叶流公式{double dfeq=0,dDotMet=0;dDotMet=ix[Order]*Velx+iy[Order]*Vely;dfeq=Den*dAlpha[Order]*(1+3*dDotMet+4.5*dDotMet*dDotMet-1.5*(Velx*Velx+Vely*Vely));return dfeq;}void prepare(float *m_df,float *u_x,float *u_y,float *dDen) //模拟前的准备工作 {for(int i=0;i<Height;i++){for(int j=0;j<Width;j++){int offset=i*Width+j;u_x[offset]=0;u_y[offset]=0;dDen[offset]=1.0;for(int l=0;l<9;l++)m_df[offset*9+l]=feq(l,Den,0,0);}}}void collide(float *bold,float *u_x,float *u_y,float *dDen)   //碰撞{for(int i=0;i<Height;i++){for(int j=0;j<Width;j++){int offset=i*Width+j;for (int l=0;l<9;l++) //碰撞 bold[offset*9+l]=bold[offset*9+l]-dRevTau*(bold[offset*9+l]-feq(l,dDen[offset],u_x[offset],u_y[offset]));for (int iPk=0;iPk<9;iPk++) //体力驱动bold[offset*9+iPk]=bold[offset*9+iPk]+Gamma[iPk]*(ix[iPk]*ForceX+iy[iPk]*ForceY);}}}void fluid_stream(float *m_df)  //流动{float *fluid_temp=(float*)malloc(Width*Height*sizeof(float)*9);for(int i=0;i<Height;i++){for(int j=0;j<Width;j++){int offset=i*Width+j;for(int l=0;l<9;l++){ int tx=remainderX(j+ix[l]),ty=remainderY(i+iy[l]);//前一个传递点的坐标int toffset=tx + ty * Width;fluid_temp[offset*9+l]=m_df[toffset*9+l];}if(i==0){fluid_temp[offset*9+1]=m_df[offset*9+3];fluid_temp[offset*9+5]=m_df[offset*9+7];fluid_temp[offset*9+8]=m_df[offset*9+6];}else if(i==Height-1){fluid_temp[offset*9+3]=m_df[offset*9+1];fluid_temp[offset*9+7]=m_df[offset*9+5];fluid_temp[offset*9+6]=m_df[offset*9+8];}}}for(int i=0;i<Height;i++){for(int j=0;j<Width;j++){int offset=i*Width+j;for(int l=0;l<9;l++){m_df[offset*9+l]=fluid_temp[offset*9+l];}}}free(fluid_temp);}void macro_calculate(float *m_df,float *u_x,float *u_y,float *den)   //计算宏观速度{for(int i=0;i<Height;i++){for(int j=0;j<Width;j++){int offset=i*Width+j;den[offset]=0;u_x[offset]=0;u_y[offset]=0;for(int l=0;l<9;l++){den[offset]+=m_df[offset*9+l];u_x[offset]+=m_df[offset*9+l]*ix[l];u_y[offset]+=m_df[offset*9+l]*iy[l];}}}}void save_data(float *m_df,float *u_x,float *u_y,float *den,int index){char namestr[30];sprintf(namestr,"data.txt");FILE *f=fopen(namestr,"w+");for(int i=0;i<Height;i++){      fprintf(f,"%d %f\n",i,u_x[Width*i+Width/2]);}fclose(f);}void make_image(unsigned char *image,float *u_x){for(int i=0;i<Height;i++){float currentValue=u_x[ i * Width + Width/2 ]; for(int j=0;j<Width * ImageScale;j++){int len=(int)(currentValue * ImageScaleU_x);int offset = i * Width * ImageScale + j;if(j>len){image[offset*4+0]=255;image[offset*4+1]=255;image[offset*4+2]=255;image[offset*4+3]=255;}else {   image[offset*4+0]=255;image[offset*4+1]=0;image[offset*4+2]=0;image[offset*4+3]=255;} }}} void anim_frame(DataBlock *d, int ticks){clock_t s=clock();collide(d->m_df,d->u_x,d->u_y,d->dDen);fluid_stream(d->m_df);macro_calculate(d->m_df,d->u_x,d->u_y,d->dDen);clock_t t=clock();cout<<ticks<<"FPS "<<t-s<<"ms"<<endl;if(ticks%100==0)save_data(d->m_df,d->u_x,d->u_y,d->dDen,ticks);    make_image(d->bitmap->get_ptr(),d->u_x);}void anim_end(DataBlock *d){free(d->m_df);free(d->u_x);free(d->u_y);free(d->dDen);}int main(){DataBlock data;CPUAnimBitmap bitmap(Width * ImageScale,Height,&data);data.bitmap=&bitmap;data.m_df=(float*)malloc(Width*Height*sizeof(float)*9);data.u_x=(float*)malloc(Width*Height*sizeof(float));data.u_y=(float*)malloc(Width*Height*sizeof(float));data.dDen=(float*)malloc(Width*Height*sizeof(float));prepare(data.m_df,data.u_x,data.u_y,data.dDen);bitmap.anim_and_exit((void(*)(void*,int))anim_frame,(void(*)(void*))anim_end);return 0;}

实验截图（在对501*501大小的方格子中进行“碰撞-流动-宏观量计算”的时间大概在160-175ms之间）：

  在做了GPU的运算以后在相同情况下，速度有所提升，但是我只是简单的运用带CUDA，感觉在效率上还会有很大的提升空间。

  下面是GPU的代码：

#include "cuda_runtime.h"#include "device_launch_parameters.h"#include "../common/book.h"#include "../common/cpu_anim.h"#include <stdio.h>#include <string.h>#define Width 501#define Height 501#define ImageScale 2          //图片宽度扩宽倍数#define ImageScaleU_x 20000   //水平速度放大倍数#define remainderX(x) ((x + Width) % Width)    //周期性边界#define remainderY(x) ((x + Height) % Height)  //周期性边界#define W0 0.4444444444444444444444444 //(4.0/9.0)#define W1 0.1111111111111111111111111 //(1.0/9.0)#define W2 0.02777777777777777777777778 //(1.0/36.0)__device__ double ForceX=1E-6;__device__ double ForceY=0;__device__ double Gamma[9] ={ 0, 1.0/3, 1.0/3, 1.0/3, 1.0/3, 1.0/12, 1.0/12, 1.0/12, 1.0/12 };__device__ double dAlpha[]={W0,W1,W1,W1,W1,W2,W2,W2,W2};__device__ int ix[9]={0, 0,1,0,-1, 1,1,-1,-1};__device__ int iy[9]={0,-1,0,1, 0,-1,1, 1,-1};__device__ double Den=1.0;__device__ double dRevTau=1.0/0.75;       //弛豫时间__device__ int inverse[9]={0,3,4,1,2,7,8,5,6};//相反方向struct DataBlock{  float *dev_boxNew;  //设备上每个点的数据  float *dev_boxOld;  //设备上每个点的数据  float *u_x;     //host上x方向上的宏观速度  float *u_y;     //host上y方向上的宏观速度  float *dDen;    //每个点的密度  CPUAnimBitmap *bitmap;};__device__ double feq(int Order,double Den,double Velx,double Vely)  //泊肃叶流公式{double dfeq=0,dDotMet=0;dDotMet=ix[Order]*Velx+iy[Order]*Vely;dfeq=Den*dAlpha[Order]*(1+3*dDotMet+4.5*dDotMet*dDotMet-1.5*(Velx*Velx+Vely*Vely));return dfeq;} __global__ void prepare(float *bold,float *u_x,float *u_y,float *den)  //模拟前的准备{int offset = threadIdx.x + blockIdx.x * blockDim.x;    for(int i=0;i<9;i++){ bold[offset*9+i]=feq(i,Den,0,0);}den[offset]=1.0;u_x[offset]=0;u_y[offset]=0;}__global__ void collide(float *bold,float *u_x,float *u_y,float *dDen)   //碰撞{   int offset = threadIdx.x + blockIdx.x * blockDim.x;for (int i=0;i<9;i++)    bold[offset*9+i]=bold[offset*9+i]-dRevTau*(bold[offset*9+i]-feq(i,dDen[offset],u_x[offset],u_y[offset]));for (int iPk=0;iPk<9;iPk++) //体力驱动bold[offset*9+iPk]=bold[offset*9+iPk]+Gamma[iPk]*(ix[iPk]*ForceX+iy[iPk]*ForceY);}__global__ void fluid_stream(float *bnew,float *bold)  //流动{   int offset = threadIdx.x + blockIdx.x * blockDim.x; for(int i=0;i<9;i++){  int tx=remainderX(threadIdx.x + ix[i]),ty=remainderY(blockIdx.x + iy[i]);//前一个传递点的坐标int toffset=tx + ty *  blockDim.x;bnew[offset*9+i]=bold[toffset*9+i];}//边界条件、半程反弹if(blockIdx.x==0){bnew[offset*9+1]=bold[offset*9+3];bnew[offset*9+5]=bold[offset*9+7];bnew[offset*9+8]=bold[offset*9+6];}if(blockIdx.x==Height-1){bnew[offset*9+3]=bold[offset*9+1];bnew[offset*9+7]=bold[offset*9+5];bnew[offset*9+6]=bold[offset*9+8];}}__global__ void macro_calculate(float *box,float *u_x,float *u_y,float *den)   //计算宏观速度{   int offset = threadIdx.x + blockIdx.x * blockDim.x;u_x[offset]=0;u_y[offset]=0;den[offset]=0;for(int i=0;i<9;i++){  u_x[offset]+=box[offset*9+i]*ix[i];u_y[offset]+=box[offset*9+i]*iy[i];den[offset]+=box[offset*9+i];}}void make_image(unsigned char *image,float *u_x){for(int i=0;i<Height;i++){float currentValue=u_x[ i * Width + Width/2 ]; for(int j=0;j<Width * ImageScale;j++){int len=(int)(currentValue * ImageScaleU_x);int offset = i * Width*2 + j;if(j>len){image[offset*4+0]=255;image[offset*4+1]=255;image[offset*4+2]=255;image[offset*4+3]=255;}else {   image[offset*4+0]=255;image[offset*4+1]=0;image[offset*4+2]=0;image[offset*4+3]=255;} }}}void save_data(DataBlock *dev,int ticks){float *m=(float*)malloc(dev->bitmap->image_size());HANDLE_ERROR(cudaMemcpy(m,dev->u_x,dev->bitmap->image_size(),cudaMemcpyDeviceToHost));   char namestr[30];sprintf(namestr,"temp/temp_%d.txt",ticks);   FILE *f=fopen(namestr,"w+");   for(int i=0;i<Height;i++)   {    fprintf(f,"%f %f\n",m[i*Width],m[i*Width + Width/2]); //if( i%Width==0 && i>0 )  fprintf(f,"\n");   }   fclose(f);   free(m);}void cleanup(DataBlock* d){cudaFree(d->dev_boxNew);    cudaFree(d->dev_boxOld);cudaFree(d->u_x);cudaFree(d->u_y);cudaFree(d->dDen);}void generate_frame(DataBlock *d,int ticks)        { clock_t s=clock();   collide<<<Height,Width>>>(d->dev_boxOld,d->u_x,d->u_y,d->dDen);   //碰撞   fluid_stream<<<Height,Width>>>(d->dev_boxNew,d->dev_boxOld);      //流动   macro_calculate<<<Height,Width>>>(d->dev_boxNew,d->u_x,d->u_y,d->dDen); //计算宏观量   float *t=d->dev_boxNew;d->dev_boxNew=d->dev_boxOld;d->dev_boxOld=t;   //交换新旧状态      clock_t c=clock();   printf("%d FPS %dms ",ticks,c-s);   if( ticks%1000==0 )  save_data(d,0);      float *temp_x=(float*)malloc(Width*Height*ImageScale*sizeof(float));      cudaMemcpy(temp_x,d->u_x,d->bitmap->image_size(),cudaMemcpyDeviceToHost);   clock_t nn=clock();   printf(" %dms ",nn-c);   make_image(d->bitmap->get_ptr(),temp_x);  //生成图片      free(temp_x);      clock_t m=clock();   printf(" %dms\n",m-nn);}int main(){    DataBlock data;CPUAnimBitmap bitmap(Width * ImageScale ,Height,&data);data.bitmap=&bitmap;HANDLE_ERROR(cudaMalloc((void**)&data.dev_boxNew,bitmap.image_size() * 9));HANDLE_ERROR(cudaMalloc((void**)&data.dev_boxOld,bitmap.image_size() * 9));HANDLE_ERROR(cudaMalloc((void**)&data.u_x,bitmap.image_size() ));HANDLE_ERROR(cudaMalloc((void**)&data.u_y,bitmap.image_size() ));HANDLE_ERROR(cudaMalloc((void**)&data.dDen,bitmap.image_size() ));prepare<<<Height,Width>>>(data.dev_boxOld,data.u_x,data.u_y,data.dDen);     save_data(&data,0);    bitmap.anim_and_exit((void(*)(void*,int))generate_frame,(void(*)(void*))cleanup);    return 0;}

下面是实验的截图（cuda在“碰撞-流动-宏观量计算”的执行时间41ms，执行时间是在CPU上执行的1/4这样）：

这个效率其实不是很高。

CPU项目文件下载地址：http://download.csdn.net/detail/tanpan004/6724617

CUDA项目文件下载地址：http://download.csdn.net/detail/tanpan004/6724621