Why im getting 1.#qnan

Question

mdgowhar 0 Newbie Poster

11 Years Ago

#include<iostream>
#include<fstream>
#include<ctime>
#include<cmath>


   using namespace std;

     clock_t start, end;
     double cpu_time_used;

int main()
{
      int m,i,j,n,k;
      double P[150][150], tou1, cx[8],cy[8], U[150][150],V[150][150];
      double u,x[150],y[150],dt=1.0, T[150][150],S[150][150],t=0.0;

cout<<"enter the units in X x Y"<<endl;
cin>>m>>n;



class pdf
        {
        public: float feq[150][150], f[150][150];
        };
pdf f0, f1,f2,f3,f4,f5,f6,f7,f8;

// For X
for (i=0;i<=m;i++)
{
    x[i]= i*1.0/(m-1);
    cout<<x[i]<<endl;
}
//For Y
for (i=n;i>=0;i--)
{

    y[i]=i*1.0/(n-1);
    cout<<y[i]<<endl;
}


//Velocity of Lattice



double mui;
cout<<"Enter Mui value"<<endl;
cin>>mui;
//Tou
tou1=((6.0*mui)+1.0)*0.5;

if (tou1>=0.55)
{
cout<<"Its Higher than Critical Value"<<endl;
}
else 
cout<<"Restart the Program & change the Mui"<<endl;


//Velocity of LID
cout<<"Enter the U velocity"<<endl;
cin>>u;

//Velocity of Lattice
cx[0]=0.0;  cy[0]=0.0;
cx[1]=1.0;  cy[1]=0.0;
cx[2]=0.0;  cy[2]=1.0;
cx[3]=-1.0; cy[3]=0.0;
cx[4]=0.0;  cy[4]=-1.0;
cx[5]=1.0;  cy[5]=1.0;
cx[6]=-1.0; cy[6]=1.0;
cx[7]=-1.0; cy[7]=-1.0;
cx[8]=1.0;  cy[8]=-1.0;

double Re;
Re=(u*m)/mui;
cout<<" Reynolds Numb:"<<Re<<endl;


for ( j=0; j<=n; j++)
{
for ( i=0 ; i<=m; i++)
{
    P[i][j]=1.0;
    U[i][j]=0.0;
    V[i][j]=0.0;
    T[i][j]=0.0;
    S[i][j]=0.0;
}}



//Weight
double w[10];
 w[0]=4.0/9.0;
       cout<<"Weight intiallized"<<endl;   
      for(i=1;i<=4;i++)
      {
                       w[i]=1.0/9.0;
                       w[i+4]=1.0/36.0;

      }

  for(i=0;i<=m;i++)
      {
U[i][n-1]=u;
}



      //initial condition
 for(j=0;j<=n;j++)
      {
      for(i=0;i<=m;i++)
      {

      f0.feq[i][j]=(w[0]*P[i][j])*(1.0-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f1.feq[i][j]=(w[1]*P[i][j])*(1.0+(3.0*U[i][j])+(4.5*(U[i][j]*U[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f2.feq[i][j]=(w[2]*P[i][j])*(1.0+(3.0*V[i][j])+(4.5*(V[i][j]*V[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f3.feq[i][j]=(w[3]*P[i][j])*(1.0-(3.0*U[i][j])+(4.5*(U[i][j]*U[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f4.feq[i][j]=(w[4]*P[i][j])*(1.0-(3.0*V[i][j])+(4.5*(V[i][j]*V[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f5.feq[i][j]=(w[5]*P[i][j])*(1.0+(3.0*(U[i][j]+V[i][j]))+(4.5*((U[i][j]+V[i][j])*(U[i][j]+V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f6.feq[i][j]=(w[6]*P[i][j])*(1.0+(3.0*(-U[i][j]+V[i][j]))+(4.5*((-U[i][j]+V[i][j])*(-U[i][j]+V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f7.feq[i][j]=(w[7]*P[i][j])*(1.0+(3.0*(-U[i][j]-V[i][j]))+(4.5*((U[i][j]+V[i][j])*(U[i][j]+V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f8.feq[i][j]=(w[8]*P[i][j])*(1.0+(3.0*(U[i][j]-V[i][j]))+(4.5*((U[i][j]-V[i][j])*(U[i][j]-V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      }
      }

 for(j=0;j<=n;j++)
      {
      for(i=0;i<=m;i++)
      { 
      f0.f[i][j]=0.0;
      f1.f[i][j]=0.0;
      f2.f[i][j]=0.0;
      f3.f[i][j]=0.0;
      f4.f[i][j]=0.0;
      f5.f[i][j]=0.0;
      f6.f[i][j]=0.0;
      f7.f[i][j]=0.0;
      f8.f[i][j]=0.0;
      } 
      }




int it,ctr=0;

//For Iteration
cout<<"Enter the Number of Iterations"<<endl;
cin>>it;
start = clock();

//*********************************************************************************************

    //time marching steps
for (int z=0; z<=it;z++)


{
      cout<<"Governing"<<endl;  
      for(i=0;i<m;i++)
      { for(j=0;j<n;j++)
      { 
            f0.f[i][j] = (1 - dt/tou1)*f0.f[i][j] + (dt/tou1)*f0.feq[i][j];
            f1.f[i][j] = (1 - dt/tou1)*f1.f[i][j] + (dt/tou1)*f1.feq[i][j];
            f2.f[i][j] = (1 - dt/tou1)*f2.f[i][j] + (dt/tou1)*f2.feq[i][j];
            f3.f[i][j] = (1 - dt/tou1)*f3.f[i][j] + (dt/tou1)*f3.feq[i][j];
            f4.f[i][j] = (1 - dt/tou1)*f4.f[i][j] + (dt/tou1)*f4.feq[i][j];
            f5.f[i][j] = (1 - dt/tou1)*f5.f[i][j] + (dt/tou1)*f5.feq[i][j];
            f6.f[i][j] = (1 - dt/tou1)*f6.f[i][j] + (dt/tou1)*f6.feq[i][j];
            f7.f[i][j] = (1 - dt/tou1)*f7.f[i][j] + (dt/tou1)*f7.feq[i][j];
            f8.f[i][j] = (1 - dt/tou1)*f8.f[i][j] + (dt/tou1)*f8.feq[i][j];              
                         } }



    //streaming

   cout<<"streaming"<<endl;
  for(j=0;j<n;j++)
            for(i=n-1;i>0;i--)
            f1.f[i][j] = f1.f[i-1][j];

        for(i=0;i<n;i++)
            for(j=n-1;j>0;j--)
            f2.f[i][j] = f2.f[i][j-1];


        for(i=n-1;i>0;i--)
            for(j=n-1;j>0;j--)
            f5.f[i][j] = f5.f[i-1][j-1];

        for(i=0;i<n-1;i++)
            for(j=n-1;j>0;j--)
            f6.f[i][j] = f6.f[i+1][j-1];

        for(i=n-1;i>0;i--)
            for(j=0;j<n-1;j++)
            f8.f[i][j] = f8.f[i-1][j+1];  

        for(i=0;i<n;i++)
            for(j=0;j<n;j++)
        {
            f0.f[i][j] = f0.f[i][j];
            if(i!=(n-1))
            f3.f[i][j] = f3.f[i+1][j];   
            if(j!=(n-1))
            f4.f[i][j] = f4.f[i][j+1];    
            if(i!=(n-1) && j!=(n-1))
            f7.f[i][j] = f7.f[i+1][j+1]; 
        }       



cout<<"boundary"<<endl;
      //boundary conditions
    for(j=0;j<n;j++)
     {
     //west                
     f1.f[0][j]=f3.f[0][j];
     f5.f[0][j]=f7.f[0][j];
     f8.f[0][j]=f6.f[0][j];

     //east
     f3.f[m-1][j]=f1.f[m-1][j];
     f6.f[m-1][j]=f8.f[m-1][j];
     f7.f[m-1][j]=f5.f[m-1][j];
     }

     for(i=0;i<m;i++)
     {
     //south
     f2.f[i][0]=f4.f[i][0];
     f5.f[i][0]=f7.f[i][0];
     f6.f[i][0]=f8.f[i][0];
     }

      for(i=1;i<m-1;i++)
      {
     //north

    t=(f0.f[i][n-1]+f1.f[i][n-1]+f3.f[i][n-1])+2*(f2.f[i][n-1]+f5.f[i][n-1]+f6.f[i][n-1]);
     f4.f[i][n-1]=f2.f[i][n-1];
     f7.f[i][n-1]=f5.f[i][n-1]-(t*u/6.0);
     f8.f[i][n-1]=f6.f[i][n-1]+(t*u/6.0);
     }



     for(j=0;j<=n;j++)
     { for(i=0;i<=m;i++)
     {

     P[i][j]=f0.f[i][j]+f1.f[i][j]+f2.f[i][j]+f3.f[i][j]+f4.f[i][j]+f5.f[i][j]+f6.f[i][j]+f7.f[i][j]+f8.f[i][j];
      }
     }
 //New Velocity
 for (i=1; i<m;i++)
P[i][n-1]=(f0.f[i][n-1]+f1.f[i][n-1]+f3.f[i][n-1])+2*(f2.f[i][n-1]+f5.f[i][n-1]+f6.f[i][n-1]);

for(j=1;j<n;j++)
{
for (i=1; i<m;i++)
{
T[i][j]=(f0.f[i][j]*cx[0])+(f1.f[i][j]*cx[1])+(f2.f[i][j]*cx[2])+(f3.f[i][j]*cx[3])+(f4.f[i][j]*cx[4])+(f5.f[i][j]*cx[5])+(f6.f[i][j]*cx[6])+(f7.f[i][j]*cx[7])+(f8.f[i][j]*cx[8]);
S[i][j]=(f0.f[i][j]*cy[0])+(f1.f[i][j]*cy[1])+(f2.f[i][j]*cy[2])+(f3.f[i][j]*cy[3])+(f4.f[i][j]*cy[4])+(f5.f[i][j]*cy[5])+(f6.f[i][j]*cy[6])+(f7.f[i][j]*cy[7])+(f8.f[i][j]*cy[8]);
U[i][j]=T[i][j]/P[i][j];
V[i][j]=S[i][j]/P[i][j];
}}

cout<<"new Equ"<<endl;
//calculation of new feq values
        for(j=0;j<=n;j++)
      {
      for(i=0;i<=m;i++)
      {
      f0.feq[i][j]=(w[0]*P[i][j])*(1.0-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f1.feq[i][j]=(w[1]*P[i][j])*(1.0+(3.0*U[i][j])+(4.5*(U[i][j]*U[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f2.feq[i][j]=(w[2]*P[i][j])*(1.0+(3.0*V[i][j])+(4.5*(V[i][j]*V[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f3.feq[i][j]=(w[3]*P[i][j])*(1.0-(3.0*U[i][j])+(4.5*(U[i][j]*U[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f4.feq[i][j]=(w[4]*P[i][j])*(1.0-(3.0*V[i][j])+(4.5*(V[i][j]*V[i][j]))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f5.feq[i][j]=(w[5]*P[i][j])*(1.0+(3.0*(U[i][j]+V[i][j]))+(4.5*((U[i][j]+V[i][j])*(U[i][j]+V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f6.feq[i][j]=(w[6]*P[i][j])*(1.0+(3.0*(-U[i][j]+V[i][j]))+(4.5*((-U[i][j]+V[i][j])*(-U[i][j]+V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f7.feq[i][j]=(w[7]*P[i][j])*(1.0+(3.0*(-U[i][j]-V[i][j]))+(4.5*((U[i][j]+V[i][j])*(U[i][j]+V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      f8.feq[i][j]=(w[8]*P[i][j])*(1.0+(3.0*(U[i][j]-V[i][j]))+(4.5*((U[i][j]-V[i][j])*(U[i][j]-V[i][j])))-(1.5*((U[i][j]*U[i][j])+(V[i][j]*V[i][j]))));
      }
      }


    ctr++; 
cout<<"Iternation No."<<ctr<<endl;     
}//time stepping ends



ofstream B,A;
A.open("2D_Uniform_LBM.plt");
      A<<"TITLE=\"2D\""<<endl;
      A<<"VARIABLES=\"X\",\"Y\",\"T\""<<endl;
      A<<"ZONE T=\"BLOCK\",I="<<n<<",J="<<m<<",F=POINT"<<endl;
B.open("2D_Uniform_LBM.txt");





for (j=n-1; j>=0; j--)
{
for (i=0; i<m;i++)
{
    cout<<U[i][j]<<"  ";
    B<<U[i][j]<<"  ";
     A<<x[i]<<"\t"<<y[j]<<"\t"<<U[i][j]<<endl;

}cout<<endl;
B<<endl;}
      end = clock();
     cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

     cout<<"The CPU time is "<< cpu_time_used<<endl;
     B<<"The CPU time is "<< cpu_time_used<<endl;
     B<<"Iterations:  "<<ctr<<endl;
B.close();
A.close();

system("pause");
}

c++

2 Contributors
3 Replies
157 Views
2 Hours Discussion Span
Latest Post 11 Years Ago Latest Post by mdgowhar

Moschops 683 Practically a Master Poster

11 Years Ago

cx[8]=1.0; cy[8]=-1.0;
cx[8] does not exist. It stops at 7. Same for cy[8].

You're doing too many iterations. You need to stop your loop earlier.

Reply to this topic

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.

mdgowhar 0 Newbie Poster · Answer 1 · 2013-09-22T09:13:50+00:00

take mui=0.01;
u=0.1
m=n=50
iteration = 1000

mdgowhar 0 Newbie Poster · Answer 2 · 2013-09-22T11:39:18+00:00

mdgowhar 0 Newbie Poster

11 Years Ago

thanks Moschops!!
:)