急问下，只能在这发！！
/******* LBM for Driven Cavity Flow with Two Lids *******/#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define M 128 // grid number in y-direction
#define N 128 // grid number in x-direction
#define M1 (M+1) // grid number in y-direction including two ends
#define N1 (N+1) // grid number in x-direction including two ends
#define Lx 1.0 // length of cavity in x-diredtion
#define Ly (Lx*M/N) // length of cavity in y-diredtion
#define U0 1.0 // Velocity in x-direction in top lid
#define V0 0.0 // Velocity in y-direction in top lid
#define SU0 0.0 // Velocity in x-direction in bottom lid
#define SV0 0.0 // Velocity in y-direction in bottom lid
#define H1 0 // grid number in 2-lids with solid boundary
#define M0 50000 // Trace Data Length
#define Q 9 // Number of discrete velocity directions of D2Q9
void lbini(void);
void Evol(void);
void exbound(void);
float feq(int k, float U, float V, float P);
void datadeal(void);
void data_read(void);
float dx,dt,c,Eps,Lam,Mu,rc,rcc;
float w;
float f[M1][N1][Q],F[M1][N1][Q];
float u[M1][N1],v[M1][N1],p[M1][N1];
float E_tr[M0],u_tr[M0],v_tr[M0],p_tr[M0];
int e[Q][2];
float rRe,Re;
int main()
{
int i,j,k,ii,m,mmax,mm,readdata;
float r,u0,v0,err=1.0,u1,v1,err1=1.0,E;
e[0][0]=e[0][1]=0;
e[1][0]=1; e[1][1]=0;
e[2][0]=0; e[2][1]=1;
e[3][0]=-1; e[3][1]=0;
e[4][0]=0; e[4][1]=-1;
e[5][0]=1; e[5][1]=1;
e[6][0]=-1; e[6][1]=1;
e[7][0]=-1; e[7][1]=-1;
e[8][0]=1; e[8][1]=-1;
printf("input Re:\n");
scanf("%f",&Re);
printf("input w:\n");
scanf("%f",&w);
dx=Lx/N;
// rRe=1.0/Re;
rRe=1.0*U0/Re; // Re=L*U0/rRe
r=(2-w)/(w+w);
dt = r*dx*dx/(3*rRe);
c=dx/dt;
rc=1.0/c;
rcc=1.0/(c*c);
Eps=5.0/12*rcc; // eps =la+mu; 0.5=la+2*mu
Mu=0.5*rcc-Eps; // mu=0.5-eps; la=2*eps-0.5
Lam=2*Eps-0.5*rcc;
lbini();
m=0;
printf("w=%f, c=%f dt=%f\n", w,c,dt);
printf("Read Data? (yes=1 no=0)\n");
scanf("%d",&readdata);
if(readdata)
{
data_read();
printf("Re same? (yes=1 no=0)\n");
scanf("%d",&readdata);
if(!readdata)
{
for(j=0;j<M1;j++) for(i=0;i<N1;i++) for(k=1;k<Q;k++)
f[j][i][k]=feq(k,u[j][i],v[j][i],p[j][i]);
}
}
// u0=u[M-2][N/2]; v0=v[M-2][N/2];
// u1=u[M/2][N/2]; v1=v[M/2][N/2];
AA:
printf("input iteration number mmax:\n");
scanf("%d",&mmax);
printf("w= %f, c=%f dt=%f\n",w,c,dt);
u0=u[M-2][N/2]; v0=v[M-2][N/2];
u1=u[M/2][N/2]; v1=v[M/2][N/2];
ii=0;
while(err>1.0e-8&&err1>5.0e-8&&ii<mmax)
// while(ii<mmax)
{
m++;
Evol();
// t +=dt;
exbound();
ii++;
if(m%100==0)
{
err=fabs(u[M-2][N/2]-u0)+fabs(v[M-2][N/2]-v0);
u0=u[M-2][N/2];v0=v[M-2][N/2];
printf("m=%d err=%e u[M-2][N2]=%e v[M-2][N2]=%e\n", m,err,u0,v0);
}
E=0.0;
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
{
E +=u[j][i]*u[j][i]+v[j][i]*v[j][i];
}
E=sqrt(E);
E_tr[m-1]=E;
u_tr[m-1]=u[M/2][N/2];
v_tr[m-1]=v[M/2][N/2];
p_tr[m-1]=p[M/2][N/2];
if(m%2000==0)
{
printf("\n");
err1=fabs(u[M/2][N/2]-u1)+fabs(v[M/2][N/2]-v1);
u1=u[M/2][N/2];v1=v[M/2][N/2];
printf("m=%d err1=%e u[M2][N2]=%e v[M2][N2]=%e\n", m,err1,u1,v1);
printf("\n");
datadeal();
}
}
datadeal();
printf("w= %f, c=%f dt=%f\n",w,c,dt);
printf("Continue? (yes=1 no=0)\n");
scanf("%d",&mm);
if(mm) goto AA;
return 0;
}
void lbini()
{
int i,j,k;
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
{u[j][i]=v[j][i]=0.0;
p[j][i]=0.0;
}
for(i=H1;i<=N-H1;i++)
{
u[M][i]=U0; v[M][i]=V0; // velocity in top lid
u[0][i]=SU0; v[0][i]=SV0; // velocity in bottom lid
}
/* for(i=0;i<=N;i++)
{
u[M][i]=16.0*U0*i*i*dx*dx*(1.0-i*dx)*(1.0-i*dx); // velocity in top lid for regularized driven flow
}
*/
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
{
for(k=1;k<Q;k++)
f[j][i][k]=feq(k,u[j][i],v[j][i],p[j][i]);
}
}
void Evol()
{
int i,j,k,id,jd;
float A,FM,s0,tmp,f57,f68;
tmp=1.0-w;
for(j=1;j<M;j++)
{
for(i=1;i<N;i++)
{
for(k=1;k<Q;k++)
{
jd=j-e[k][1]; id =i-e[k][0];
FM=feq(k,u[jd][id],v[jd][id],p[jd][id]);
F[j][i][k]=tmp*f[jd][id][k]+w*FM;
}
}
}
tmp=0.25/Eps;
for(j=1;j<M;j++)
{
for(i=1;i<N;i++)
{
u[j][i]=v[j][i]=p[j][i]=A=0.0;
for(k=1;k<Q;k++) A+=F[j][i][k];
for(k=1;k<Q;k++) f[j][i][k]=F[j][i][k];
f57=f[j][i][5]-f[j][i][7];
f68=f[j][i][6]-f[j][i][8];
u[j][i]=f[j][i][1]-f[j][i][3]+f57-f68;
v[j][i]=f[j][i][2]-f[j][i][4]+f57+f68;
/* for(k=1;k<Q;k++)
{
u[j][i] +=(f[j][i][k]*e[k][0]);
v[j][i] +=(f[j][i][k]*e[k][1]);
}
*/
s0=-2.0/3*(u[j][i]*u[j][i]+v[j][i]*v[j][i]);
u[j][i] *= c;
v[j][i] *= c;
p[j][i] = A+s0;
p[j][i] *= tmp;
}
}
}
float feq(int k,float U, float V, float P)
{
float su,uv,eu,x;
eu=(e[k][0]*U+e[k][1]*V)*rc;
uv=U*U+V*V;
su=3*eu+4.5*eu*eu-1.5*uv*rcc;
/* if(k==0) x=-4*Eps*P+su*4.0/9;*/
if(k<5) x=Lam*P+su/9;
else x=Mu*P+su/36;
return x;
}
/*
void datadeal()
{
int i,j;
FILE *fp;
if( ( fp=fopen("u","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<=M;j++)
{
for (i=0; i<=N; i++) fprintf(fp,"%e ",u[j][i]);
fprintf(fp,"\n");
}
fclose(fp);
if( ( fp=fopen("v","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<=M;j++)
{
for (i=0; i<=N; i++) fprintf(fp,"%e ",v[j][i]);
fprintf(fp,"\n");
}
fclose(fp);
if( ( fp=fopen("p","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<=M;j++)
{
for (i=0; i<=N; i++) fprintf(fp,"%e ",p[j][i]);
fprintf(fp,"\n");
}
fclose(fp);
if(( fp=fopen("E_tr","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<M0;j++)
{
fprintf(fp,"%f ",E_tr[j]);
fprintf(fp,"\n");
}
fclose(fp);
if(( fp=fopen("u_tr","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<M0;j++)
{
fprintf(fp,"%f ",u_tr[j]);
fprintf(fp,"\n");
}
fclose(fp);
if(( fp=fopen("v_tr","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<M0;j++)
{
fprintf(fp,"%f ",v_tr[j]);
fprintf(fp,"\n");
}
fclose(fp);
if(( fp=fopen("p_tr","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<M0;j++)
{
fprintf(fp,"%f ",p_tr[j]);
fprintf(fp,"\n");
}
fclose(fp);
fp=fopen("ut","wb");
fwrite(u,sizeof(float),M1*N1,fp);
fclose(fp);
fp=fopen("vt","wb");
fwrite(v,sizeof(float),M1*N1,fp);
fclose(fp);
fp=fopen("pt","wb");
fwrite(p,sizeof(float),M1*N1,fp);
fclose(fp);
fp=fopen("ft","wb");
fwrite(f,sizeof(float),M1*N1*Q,fp);
fclose(fp);
}
*/
void datadeal(void)
{
int x,y;
float tmp;
char var_name[5][32];
FILE *fp=NULL;
strcpy(var_name[0],"X");
strcpy(var_name[1],"Y");
strcpy(var_name[2],"U");
strcpy(var_name[3],"V");
strcpy(var_name[4],"P");
fp=fopen("U.plt","wb");
TecplotHead(fp,5,var_name,"zone1",M1,N1,1);
for(x=1;x<=M1;x++)
for(y=1;y<=N1;y++)
{
tmp=y;
fwrite(&tmp,1,sizeof(float),fp);
}
for(x=1;x<=M1;x++)
{
tmp=x;
for(y=1;y<=N1;y++)
fwrite(&tmp,1,sizeof(float),fp);
}
fwrite(u,M1 * N1,sizeof(float),fp);
fwrite(v,M1 * N1,sizeof(float),fp);
fwrite(p,M1 * N1,sizeof(float),fp);
fclose(fp);
}
void data_read()
{
FILE *fp;
fp=fopen("ut","rb");
fread(u,sizeof(float),M1*N1,fp);
fclose(fp);
fp=fopen("vt","rb");
fread(v,sizeof(float),M1*N1,fp);
fclose(fp);
fp=fopen("pt","rb");
fread(p,sizeof(float),M1*N1,fp);
fclose(fp);
fp=fopen("ft","rb");
fread(f,sizeof(float),M1*N1*Q,fp);
fclose(fp);
}
void exbound()
{
int i,j,k;
float FMb,FMf;
for(j=0;j<=M;j++)
{
// u[j][0]=u[j][N]=v[j][0]=v[j][N]=0.0;
p[j][0]=p[j][1];
p[j][N]=p[j][N-1];
for(k=0;k<Q;k++)
{
FMb=feq(k,u[j][0],v[j][0],p[j][0]); //Noneq expolation
FMf=feq(k,u[j][1],v[j][1],p[j][1]);
f[j][0][k]=FMb+(f[j][1][k]-FMf);
FMb=feq(k,u[j][N],v[j][N],p[j][N]);
FMf=feq(k,u[j][N-1],v[j][N-1],p[j][N-1]);
f[j][N][k]=FMb+(f[j][N-1][k]-FMf);
}
}
for(i=0;i<H1;i++)
{
u[0][i]=u[M][i]=0.0; // solid boundary in 2-lid
u[0][N-i]=u[M][N-i]=0.0;
}
for(i=H1;i<=N-H1;i++)
{
u[M][i]=U0; v[M][i]=V0; // velocity in top lid
u[0][i]=SU0; v[0][i]=SV0; // velocity in bottom lid
}
/* for(i=0;i<=N;i++)
{
u[M][i]=16.0*U0*i*i*dx*dx*(1.0-i*dx)*(1.0-i*dx); // velocity in top lid
}
*/
for(i=0;i<=N;i++)
{
p[0][i]=p[1][i];
p[M][i]=p[M-1][i];
for(k=0;k<Q;k++)
{
FMb=feq(k,u[0][i],v[0][i],p[0][i]);
FMf=feq(k,u[1][i],v[1][i],p[1][i]);
f[0][i][k]=FMb+(f[1][i][k]-FMf);
FMb=feq(k,u[M][i],v[M][i],p[M][i]);
FMf=feq(k,u[M-1][i],v[M-1][i],p[M-1][i]);
f[M][i][k]=FMb+(f[M-1][i][k]-FMf);
}
}
// u[0][0]=u[1][1]; u[0][N]=u[1][N-1];
// u[M][0]=u[M-1][1]; u[M][N]=u[M-1][N-1];
// v[0][0]=v[1][1]; v[0][N]=v[1][N-1];
// v[M][0]=v[M-1][1]; v[M][N]=v[M-1][N-1];
p[0][0]=p[1][1]; p[0][N]=p[1][N-1];
p[M][0]=p[M-1][1]; p[M][N]=p[M-1][N-1];
for(k=1;k<Q;k++)
{
FMb=feq(k,u[0][0],v[0][0],p[0][0]);
FMf=feq(k,u[1][1],v[1][1],p[1][1]);
f[0][0][k]=FMb+(f[1][1][k]-FMf);
FMb=feq(k,u[0][N],v[0][N],p[0][N]);
FMf=feq(k,u[1][N-1],v[1][N-1],p[1][N-1]);
f[0][N][k]=FMb+(f[1][N-1][k]-FMf);
FMb=feq(k,u[M][0],v[M][0],p[M][0]);
FMf=feq(k,u[M-1][1],v[M-1][1],p[M-1][1]);
f[M][0][k]=FMb+(f[M-1][1][k]-FMf);
FMb=feq(k,u[M][N],v[M][N],p[M][N]);
FMf=feq(k,u[M-1][N-1],v[M-1][N-1],p[M-1][N-1]);
f[M][N][k]=FMb+(f[M-1][N-1][k]-FMf);
}
}
/* the end */
上面的程序TecplotHead(fp,5,var_name,"zone1",M1,N1,1)函数不能识别，是什么问题？
还有void Evol()用的什么算法？？？
