Hello friends,
I am trying to implement lattice Boltzmann method for solving 2D diffusion problem for a square cavity where top and bottom boundaries are set to neumann boundary condition and left and right boundaries are set dirichlet condition. Below is the code. The iteration is getting terminated with a “nan”. I have tried a lot to deal with it, went through palabos sites,openlb sites but all went unsolving. Please any one can rescue me. It would be a great help.
#include<stdio.h>
#define nx 100
#define ny 100
void main()
{
int i,j,k,timestep=0;
FILE *a,*b,c;
a=fopen(“grid.dat”,“w”);
b=fopen(“isotherm.dat”,“w”);
c=fopen(“thermal.dat”,“w”);
float x[101],y[101],th[101][101],g[101][101][10],th0[101][101],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};
float dx,dy,taug,e1,e2,error,alpha;
float sumt=0.0,tl=1.0,tr=0.0,lx=1.0,ly=1.0;
fprintf(a,“ZONE I=101\t J=101\n”);
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]=ilx/nx;
y[j]=jly/ny;
fprintf(a,"%f\t %f\n", x[i],y[j]);
}
}
/----------------------------initial condition----------------------------/
for(j=0;j<=ny;j++)
{
for(i=1;i<nx;i++)
{
th[i][j]=1.0;
for(k=0;k<9;k++)
{
g[i][j][k]=w_g[k]*th[i][j];
//printf(“pdf=%f\n”,g[i][j][k]);
}
}
}
for(j=0;j<=ny;j++)
{
for(k=0;k<9;k++)
{
g[0][j][k]=w_g[k]tl;
//printf(“left pdf=%f\n”,g[0][j][k]);
g[nx][j][k]=w_g[k]tr;
//printf(“right pdf=%f\n”,g[nx][j][k]);
}
}
/-------------------------------main loop---------------------------------/
error=10.0;
do
{
/-------------------------Collision step-----------------------------/
for(j=0;j<=ny;j++)
{
for(i=1;i<nx;i++)
{
for(k=0;k<9;k++)
{
geq=w_g[k]*th[i][j];
g[i][j][k]=(g[i][j][k]+((geq-g[i][j][k])/taug));
}
}
}
for(j=0;j<=ny;j++)
{
for(k=0;k<9;k++)
{
geq=w_g[k]tl;
g[0][j][k]=(g[0][j][k]+((geq-g[0][j][k])/taug));
}
}
for(j=0;j<=ny;j++)
{
for(k=0;k<9;k++)
{
geq=w_g[k]tr;
g[nx][j][k]=(g[nx][j][k]+((geq-g[nx][j][k])/taug));
}
}
/---------------------------------------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])-g[0][j][3];
g[0][j][5]=tl*(w_g[5]+w_g[7])-g[0][j][7];
g[0][j][8]=tl*(w_g[8]+w_g[6])-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];
}
/-----------------------------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];
}
}
/-----------------------------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(“th=%f\n”,th[i][j]);
}
}
/-------------------------------------error checking condition------------------------------------------/
e1=e2=0.0;
for(i=1;i<nx;i++)
{
for(j=0;j<=ny;j++)
{
e1=e1+(th[i][j]-th0[i][j]);
e2=e2+th[i][j];
th0[i][j]=th[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(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©;
}
How to deal with this NAN??