Lid driiven cavity code using LBM

Hello Everyone

I am a beginner in this field and trying simple codes in LBM.I have wriiten the code but is not working and not giving suitable results. If some one can help me.
here is the code

#include<stdio.h>
#include"crun4.h"

void main()
{
FILE *a,ifp,fp3,fp4;
a=fopen(“lid.dat”,“w”);
ifp = fopen(“grid.dat”, “w”);
fp3=fopen(“uvel.dat”,“w”);
fp4 = fopen(“vvel.dat”, “w”);
fprintf(ifp,“ZONE I=200\t J=200”);
fprintf(ifp,"\n");
dx=1.0/nx;
dy=1.0/ny;
dt=dx;
tau=0.7;
nu=((2
tau-1.0)dxdx)/(6
dt);
u1=(nu
Rex);
x[0][0]=0.0;
y[0][0]=0.0;
x[nx][j]=1.0;
y[i][ny]=1.0;

/*----------------------------grid generation----------------------------*/
	for(i=1;i<nx;i++)
	{
            for(j=1;j<ny;j++)
	  {	

	
		x[i][j]=i*dx;
		y[i][j]=j*dy;
                fprintf(ifp,"%lf\t %lf\n", x[i][j],y[i][j]);
		//printf("grids are: %f\t%f\n",x[i],y[j]);
	}
}

/*----------------------------initial condition----------------------------*/

for(i=0;i<=nx;i++)
{
	for(j=0;j<ny;j++)
		{
			rho[i][j]=1.0;
			ux[i][j]=0.0;
			uy[i][j]=0.0;
			for(k=0;k<9;k++)
				{
					f[i][j][k]=w_f[k]*rho[i][j];
					fp[i][j][k]=w_f[k]*rho[i][j];

// //printf(“initialised PDF\n”);
}
}
}
for(i=1;i<nx;i++)
{
rho[i][ny]=1.0;
ux[i][ny]=u1;
uy[i][ny]=0.0;
// F2=0.01;
// for(k=0;k<9;k++)
// {
// F1=ex[k]0.1;
// feq=w_f[k]rho[i][ny](1.0+(3.0
F1)+(4.5F1F1)-(1.5*F2));
// f[i][ny][k]=feq;
// //printf("%f\n",f[i][ny][k]);
// }
}
for(i=0;i<=nx;i++)
{

      // rhon=f[i][ny][0]+f[i][ny][1]+f[i][ny][3]+2.0*((f[i][ny][2])+(f[i][ny][6])+(f[i][ny][5]));
	F2=u1*u1;
	for(k=0;k<9;k++)
	{
		F1=ex[k]*u1;
		feq=w_f[k]*rho[i][ny]*(1.0+(3.0*F1)+(4.5*F1*F1)-(1.5*F2));			
		f[i][ny][k]=feq;

// //printf("%f\n",f[i][ny][k]);
}
/f[i][ny][4]=f[i][ny][2];
f[i][ny][8]=f[i][ny][6]+(rho[i][ny]
(ux[i][ny]/6.0));
f[i][ny][7]=f[i][ny][5]+(rho[i][ny](ux[i][ny]/6.0));/

}


/*-------------------------------main loop---------------------------------*/

error=10.0;
do
{
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
E2=ux[i][j]ux[i][j]+uy[i][j]uy[i][j];
for(k=0;k<9;k++)
{
E1=ex[k]ux[i][j]+ey[k]uy[i][j];
feq=w_f[k]rho[i][j](1.0+(3.0
E1)+(4.5
E1
E1)-(1.5
E2));
f[i][j][k]=fp[i][j][k]+((feq-fp[i][j][k])/tau);
//printf("%f\n",fp[i][j][k]);
}
}
}

	/*-----------------------------Streaming step-----------------------------*/	
	for(j=0;j<=ny;j++)
{
	for(i=nx;i>=1;i--)
	{
		fp[i][j][1]=f[i-1][j][1];

	}

	for(i=0;i<nx;i++)
	{
		fp[i][j][3]=f[i+1][j][3];
	}
}
for(j=ny;j>=1;j--)							/*At j=0,we have pdfs from Boundary conditions*/
{
	for(i=0;i<=nx;i++)
	{
		fp[i][j][2]=f[i][j-1][2];
	}
	for(i=nx;i>0;i--)
	{
		fp[i][j][5]=f[i-1][j-1][5];
	}
	for(i=0;i<nx;i++)
	{
		fp[i][j][6]=f[i+1][j-1][6];
	}
}
for(j=0;j<ny;j++)							/*At j=ny,pdfs are defined*/
{
	for(i=0;i<=nx;i++)
	{
		fp[i][j][4]=f[i][j+1][4];
	}
	for(i=0;i<nx;i++)
	{
		fp[i][j][7]=f[i+1][j+1][7];
	}
	for(i=nx;i>=1;i--)
	{
		fp[i][j][8]=f[i-1][j+1][8];
	}
}


/*---------------------------------------Boundary Condition--------------------------------------*/

/*---------------------left wall boundary condition-----------------------*/
for(j=1;j<ny;j++)
{
	f[0][j][1]=fp[0][j][3];
	f[0][j][5]=fp[0][j][7];
	f[0][j][8]=fp[0][j][6];
}
/*---------------------right wall boundary condition-----------------------*/
for(j=1;j<ny;j++)
{
	f[nx][j][3]=fp[nx][j][1];
	f[nx][j][6]=fp[nx][j][8];
	f[nx][j][7]=fp[nx][j][5];
}
/*---------------------bottom wall boundary condition-----------------------*/

for(i=0;i<=nx;i++)
{
	f[i][0][2]=fp[i][0][4];
	f[i][0][5]=fp[i][0][7];
	f[i][0][6]=fp[i][0][8];
}
 

/*-----------------------------Calculation of velocity and density-----------------------------*/
for(i=1;i<nx;i++)
	{
		for(j=1;j<ny;j++)
			{
			  sum=0.0;
			   for(k=0;k<9;k++)
			   {
			     sum=sum+f[i][j][k];
			   }
			   rho[i][j]=sum;

// rho[i][j]=f[i][j][1]+f[i][j][2]+f[i][j][3]+f[i][j][4]+f[i][j][5]+f[i][j][6]+f[i][j][7]+f[i][j][8]+f[i][j][0];
// ux[i][j]=((f[i][j][1]+f[i][j][5]+f[i][j][8])-(f[i][j][3]+f[i][j][6]+f[i][j][7]))/rho[i][j];
// uy*/[i][j]=((f[i][j][2]+f[i][j][5]+f[i][j][6])-(f[i][j][4]+f[i][j][7]+f[i][j][8]))/rho[i][j];
}
}

	for(i=1;i<nx;i++)
	{
		for(j=1;j<ny;j++)
			{
			  uxsum=0.0;
			  uysum=0.0;
			   for(k=0;k<9;k++)
			   {
			     uxsum=uxsum+(f[i][j][k]*ex[k]);
			     uysum=uysum+(f[i][j][k]*ey[k]);
			   }
			   ux[i][j]=uxsum/rho[i][j];
			   uy[i][j]=uysum/rho[i][j];

			}
	}
	
/*-------------------------------------error checking condition------------------------------------------*/

e1=e2=0.0;
for(i=1;i<nx;i++)
{
	for(j=1;j<ny;j++)
	{
		e1=e1+(sqrt((ux[i][j]-ux0[i][j])*(ux[i][j]-ux0[i][j])+(uy[i][j]-uy0[i][j])*(uy[i][j]-uy0[i][j])));
		e2=e2+(sqrt((ux[i][j]*ux[i][j])+(uy[i][j]*uy[i][j])));
		ux0[i][j]=ux[i][j];
		uy0[i][j]=uy[i][j];
	}
}
error=e1/e2;
if(timestep%1000==0)
{
	printf("%d\t%0.13f\n",timestep,error);
}
	timestep++;

// if(timestep>=100000) error=0.0;
}
while(error>=1e-8);
printf("%d\t%0.13f\n",timestep,error);
/---------------------------------------------output------------------------------------------------------/
fprintf(a,“Zone I=%d\tJ=%d \n”,nx+1,ny+1 );
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
fprintf(a,"%f\t%f\t%f\t %f\n",x[i],y[j],ux[i][j],uy[i][j]);
}
}
for(j=1;j<=ny;j++)
{
fprintf(fp3,"\n%f\t%f",y[100][j],ux[100][j]);
}
for(i=1;i<=nx;i++)
{
fprintf(fp4,"\n%f\t%f",x[i][100],uy[i][100]);
}
fclose(ifp);
fclose(a);
fclose(fp3);
fclose(fp4);
//getch();
}