I have written a 2D code to solve conduction equation where the left and right walls are set to dirichlet boundary condition and bottom and top walls are set to neumann boundary condition. I am posting my code.Please help where I am going wrong because I am not getting any output.
#include<stdio.h>
#include<math.h>
#define nx 100
#define ny 100
void main()
{
int i,j,k,itrstep,timestep=0;
FILE *a,*b,c;
a=fopen(“grid.dat”,“w”);
b=fopen(“isotherm.dat”,“w”);
c=fopen(“thermal.dat”,“w”);
double x[101],y[101],th[101][101],g[101][101][10],geq;
double w_g[9]={(4.0/9.0),(1.0/9.0),(1.0/9.0),(1.0/9.0),(1.0/9.0),(1.0/36.0),(1.0/36.0),(1.0/36.0),(1.0/36.0)};
int ex[9]= {0,1,0,-1,0,1,-1,-1,1};
int ey[9]= {0,0,1,0,-1,1,1,-1,-1};
double dx,dy,E1,E2,taug,nu,F1,F2,e1,e2,error,alpha;
double sumt=0.0,tl=1.0,tr=0.0,lx=1.0,ly=1.0;
dx=lx/nx;
dy=ly/ny;
alpha=0.16;
taug=(3alpha+0.5);
/*----------------------------grid generation----------------------------*/
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
x[i]=x[i-1]*lx/nx;
y[j]=y[j-1]*ly/ny;
fprintf(a,"%f\t %f\n", x[i],y[j]);
}
}
/*----------------------------initial condition----------------------------*/
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
th[i][j]=0.0;
for(k=0;k<9;k++)
{
g[i][j][k]=w_g[k]*th[i][j];
if(i==0)
{
g[i][j][k]=w_g[k]*tl;
}
if(i==nx)
{
g[i][j][k]=w_g[k]*tr;
}
}
}
}
/*-------------------------------main loop---------------------------------*/
itrstep=5000;
do
{
/-------------------------Collision step-----------------------------/
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
for(k=0;k<9;k++)
{
geq=g[i][j][k]=w_g[k]*th[i][j];
g[i][j][k]=(g[i][j][k]+((geq-g[i][j][k])/taug));
}
}
}
/*-----------------------------Streaming step-----------------------------*/
for(j=0;j<=ny;j++)
{
for(i=nx;i>=1;i--)
{
g[i][j][1]=g[i-1][j][1];
}
for(i=0;i<nx;i++)
{
g[i][j][3]=g[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++)
{
g[i][j][2]=g[i][j-1][2];
}
for(i=nx;i>0;i--)
{
g[i][j][5]=g[i-1][j-1][5];
}
for(i=0;i<nx;i++)
{
g[i][j][6]=g[i+1][j-1][6];
}
}
for(j=0;j<ny;j++) /*At j=ny,pdfs are defined*/
{
for(i=0;i<=nx;i++)
{
g[i][j][4]=g[i][j+1][4];
}
for(i=0;i<nx;i++)
{
g[i][j][7]=g[i+1][j+1][7];
}
for(i=nx;i>=1;i--)
{
g[i][j][8]=g[i-1][j+1][8];
}
}
/*---------------------------------------Boundary Condition for g--------------------------------------*/
/*---------------------left and right wall boundary condition-----------------------*/
for(j=0;j<=ny;j++)
{
g[0][j][1]=tl*w_g[1]+w_g[3]*tl-g[0][j][3];
g[0][j][5]=tl*w_g[5]+w_g[7]*tl-g[0][j][7];
g[0][j][8]=tl*w_g[8]+w_g[6]*tl-g[0][j][6];
g[nx][j][6]=-g[nx][j][8];
g[nx][j][3]=-g[nx][j][1];
g[nx][j][7]=-g[nx][j][5];
}
/*---------------------Top and bottom adiabatic wall boundary condition----------------*/
for(i=0;i<=nx;i++)
{
g[i][ny][8]=g[i][ny-1][8];
g[i][ny][7]=g[i][ny-1][7];
g[i][ny][6]=g[i][ny-1][6];
g[i][ny][5]=g[i][ny-1][5];
g[i][ny][4]=g[i][ny-1][4];
g[i][ny][3]=g[i][ny-1][3];
g[i][ny][2]=g[i][ny-1][2];
g[i][ny][1]=g[i][ny-1][1];
g[i][ny][0]=g[i][ny-1][0];
g[i][0][8]=g[i][1][8];
g[i][0][7]=g[i][1][7];
g[i][0][6]=g[i][1][6];
g[i][0][5]=g[i][1][5];
g[i][0][4]=g[i][1][4];
g[i][0][3]=g[i][1][3];
g[i][0][2]=g[i][1][2];
g[i][0][1]=g[i][1][1];
g[i][0][0]=g[i][1][0];
}
/*-----------------------------Calculation of temperature-----------------------------*/
sumt=0.0;
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
for(k=0;k<9;k++)
{
sumt=sumt+g[i][j][k];
}
th[i][j]=sumt;
}
}
printf("%d\n",timestep);
timestep++;
}while(timestep<=itrstep);
/*---------------------------------------------output------------------------------------------------------*/
fprintf(b,"I=%d\tJ=%d \n",nx+1,ny+1);
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
fprintf(b,"%f\t%f\t%f\n",x[i],y[j],th[i][j]);
}
}
for(i=0;i<=nx;i++)
{
fprintf(c,"\n%f\t%f",x[i],th[i][ny/2]);
}
fclose(a);
fclose(b);
fclose©;
}
Earliest help is highly encouraged.Thanx.
???