2D diffusion problem using Lattice Boltzmann method

shikha5

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??

Chrisir

Is this C?

Format your code

quark

This is C

koogs

edit post, highlight code, press ctrl-o

Kindly copy it to editpad pro and then proceed please.

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.

koogs

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

quark

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

shikha5

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.

quark

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 this

/** ------------------ grid generation ------------------- */
for (i=0; i<=nx; i++)
{
  x[i] = i * lx / nx;
}

for (i=0; i<=ny; i++)
{
  y[i] = i * ly / ny;
}

I 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.

shikha5

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.

quark

Also, I have started the loop from 0 to nx so x[0],x[nx]and y[0],y[ny] have got their values.

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.

So, I would also like to ask you that my grid is generated finely as I have plotted it.

Not sure what you mean by generated finely as I have plotted it.

koogs

i tried porting this to processing

i get a NaN around

while (error>=1e-8);

error is being set by error=e1 / e2; which are both Infinity according to the debugger.