#include<iostream>
#include<cmath>
#include<cstdlib>
#include<iomanip>
#include<fstream>
#include<sstream>
#include<string>
using namespace std;
const int Q=9;
const int NX=256;
const int NY=256;
const double U=0.1;
int e[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
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],f[NX+1][NY+1][Q],F[NX
+1][NY+1][Q];
int i,j,k,ip,jp,n;
double c,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error;
void init();
double feq(int k,double rho,double u[2]);
void evolution();
void output(int m);
void Error();
int main()
{
using namespace std;
init();
for(n=0;;n++)
{
evolution();
if(n%100==0)
{
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)
{
评论0