基于 格子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