C++编程：利用格子玻尔兹曼方法模拟自然对流

# include "stdafx.h"
# include <iostream> 
# include <cmath> 
# include <cstdlib> 
# include <iomanip> 
# include <fstream> 
# include <sstream> 
# include <string> 

using namespace std;                    
const int Q=9;         
const int NX=79;       //x方向格子数
const int NY=79;         //y方向格子数
const double T2=1,T1=0;      //两端温度,这里温度定义为内能密度，即内能除以物质量


double e[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};             //直角方向 9个移动方向
double w[Q]={4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36};              //不同方向的值
double rho[NX+1][NY+1],u[NX+1][NY+1][2],u0[NX+1][NY+1][2],T[NX+1][NY+1],g[NX+1][NY+1][Q],G[NX+1][NY+1][Q]; //密度,速度,分布函数,温度场,内能分布函数                                             //两个系数和加速度
int i,j,k,ip,jp,n;
double c,Ra,Pr,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,tau_g,niu,rsk,beta,tref,temp,error;

void init();                                                                     //函数声明 
double feqf(int k, double rho, double u[2]);
double feqt(int k, double T, double u[2]);
void evolution();
void outputf(int m);
void outputt(int m);
void Error();

int main()
{
using namespace std;                                                //变量名声明
init();                                                             //初始化
for(n=0; ;n++)                                                      //循环
{
evolution();                                                   //1.演化计算
if(n%200==0)                                                   //2.每200步输出一次
{
Error();                                               //计算残值
cout<<"The"<<n<<"th computation result:"<<endl<<"The u,v of point (NX/2,NY/2) is :"<<setprecision(6)<<u[NX/2][NY/2][0]<<","<<u[NX/2][NY/2][1]<<endl;
cout<<"The max relative error of uv is:"<<setiosflags(ios::scientific)<<error<<endl;                                                           //输出中心值和残值
if(n>=1000)
{
if(n%1000==0)                                  //每1000步写出一次
outputf(n);                                //收敛判断
    outputt(n);
if(error<1.0e-6)
break;
}
}
}
  return 0;
}

void init()                                                //初始化函数
{
dx=1.0;
dy=1.0;
Lx=dx*double(NX);;                    //物理参数设置，每格长度为1.0