We are about to switch to a new forum software. Until then we have removed the registration on this forum.

- All Categories 25.7K
- Announcements & Guidelines 13
- Common Questions 30
- Using Processing 22.1K
- Programming Questions 12.2K
- Questions about Code 6.4K
- How To... 4.2K
- Hello Processing 72
- GLSL / Shaders 292
- Library Questions 4K
- Hardware, Integration & Other Languages 2.7K
- Kinect 668
- Arduino 1K
- Raspberry PI 188
- Questions about Modes 2K
- Android Mode 1.3K
- JavaScript Mode 413
- Python Mode 205
- Questions about Tools 100
- Espanol 5
- Developing Processing 548
- Create & Announce Libraries 211
- Create & Announce Modes 19
- Create & Announce Tools 29
- Summer of Code 2018 93
- Rails Girls Summer of Code 2017 3
- Summer of Code 2017 49
- Summer of Code 2016 4
- Summer of Code 2015 40
- Summer of Code 2014 22
- p5.js 1.6K
- p5.js Programming Questions 947
- p5.js Library Questions 315
- p5.js Development Questions 31
- General 1.4K
- Events & Opportunities 289
- General Discussion 365

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=(3*alpha+0.5);
/** ----------------------------grid generation--------------------- */
for (j=0; j<=ny; j++)
{
for (i=0; i<=nx; i++)
{
x[i]=i*lx/nx;
y[j]=j*ly/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 -------------- */
/** ---------------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(c);
}
```

How to deal with this NAN??

## Answers

Is this C?

Format your code

This is C

edit post, highlight code, press ctrl-o

this won't help - all the pairs of * symbols have been converted to italics markup. the only thing that will help is you fixing the post.

but given that it's C and not processing3 then we probably won't be able to help.

I have formatted the code an modified the comments so that the code is legible.

Can you give me your mail Chrisir and Quark so that I can send you the formatted code directly. Actually I am unable to format my code properly here.

`float x[101], y[101], th[101][101], g[101][101][10], th0[101][101], geq;`

In C when you declare a variable it reserves memory for that variable but does NOT initialise the variable to a valid value.

So if we consider the

`x`

array as an example it will reserve space for 101 floating point values, that's 404 bytes but the contents of this memory could be anything. Some of the array elements might actually contain valid floating point numbers but others won't so will be reported as NaN.Therefore in C it is important to make sure that all elements are initialised to some valid value for the data type.

In the grid generation you are not initializing ALL the elements in the arrays. For instance

`x[0]`

and`x[nx]`

are not being initialised. BTW you are initialising the x array ny times :( you don't need the double loop because it is not a 2D array. Try thisI don't know if this is the cause of your NaN but since your code is C and this forum is for Processing (Java) no one is likely to fix the code rather just offer suggestions.

Thanx a lot Quark or your kind suggestion. But I am also new to the coding world. So, I would also like to ask you that my grid is generated finely as I have plotted it. Also, I have started the loop from 0 to nx so x[0],x[nx]and y[0],y[ny] have got their values. So, where I am going wrong. Please enlighten me.

Sorry about that I was looking at the wrong line although my comment about the double loop is still valid.

Also my comments about making sure all variables and array elements are initialised is still true for C. BTW I just remembered C does not do array bounds checking so you need to make sure you don't iterate over non-existent array elements as well.

Do you know where the NaN is first detected?

This type of error is extremely hard to find because it is a runtime error. I suggest that you step through your code using your IDE and examine the variables to see where the problem is.

Not sure what you mean by

generated finely as I have plotted it.i tried porting this to processing

i get a NaN around

error is being set by

`error=e1 / e2;`

which are both Infinity according to the debugger.