Linux下c程序运行结果出现nan

长刺的耗子

问题：
相同的c程序，在Windows下vs2010环境下运行结果正确，输出的double型

但是，在Linux环境下，g++编译，运行结果出现nan

MMMIX

回复 1# 长刺的耗子


    代码呢？另外，C 代码你干吗用 g++ 编译？

不过，出现这个情况，基本都是由于你的代码写的有问题。

长刺的耗子

因为程序后缀是cpp，所以用的是g++，不过改成.c文件，用gcc编译也是同样的错误


#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"



//#include "stdafx.h"
//#include "math.h"
//#include "stdio.h"

#define MaxGrid 200
#define PI 3.1415926
double p[MaxGrid+2][MaxGrid+2][MaxGrid+2],u[MaxGrid+2][MaxGrid+2][MaxGrid+2],newq[MaxGrid+2][MaxGrid+2][MaxGrid+2],newc[MaxGrid+2][MaxGrid+2][MaxGrid+2];
double c[MaxGrid+2][MaxGrid+2][MaxGrid+2],fcc[MaxGrid+2][MaxGrid+2][MaxGrid+2],fc[MaxGrid+2][MaxGrid+2][MaxGrid+2],cl[MaxGrid+2][MaxGrid+2][MaxGrid+2],cs[MaxGrid+2][MaxGrid+2][MaxGrid+2];
double f[MaxGrid+2][MaxGrid+2][MaxGrid+2];
//double gp1[MaxGrid+5][MaxGrid+5],hp1[MaxGrid+5][MaxGrid+5];
int num;
int i,j,k;
double t;
///////////////////变量赋初值及无量刚化/////////////////

double deltx=0.6e-8,delty=0.6e-8,deltz=0.6e-8;
double delt=1.5e-9;

///////////////////////////////////////////////////////////////////////

double jmn=0.093;         //界面能

double R=8.31441;          //气体常数

double vm=10.547e-6;          //摩尔体积
double ke=0.14;            //平衡分配系数

double c0=2.0/100.0;
double Ds=3.0e-13;             //固相扩散率
double Dl=3.0e-9;         //液相中的溶质扩散系数
double Tm=933.3;
double T=885;//熔点温度
double x=4.5*deltx/2.0;    //界面厚度

double ka=4;                //对称性系数
double v=0.063;             //各向异性强度
double M;//=0.1561227;                  //相场参数

double cle=(933.3-T)/620;     //平衡溶质分配系数
double cse=ke*cle;
double eip0=sqrt(6*jmn*x/2.2);               //eip
double eip,eip1,eip2;
double eipx,eipx1,eipx2;
double eipy,eipy1,eipy2;
double eipz,eipz1,eipz2;


double W=6.6*jmn/x;

double jifen;
double hp;
//double hp1;     //hp的一阶导数
double gp;
///double gp1;     //gp的一阶导数
double pa;

//double fc;
//double fcc;
double fclcl;
double fcscs;
double Dp;        //  扩散系数     
double dDp;           //Dp的梯度
double dhp;           //hp的梯度
double dcl;           //cl的梯度
double dcs;           //cs的梯度
double dlncl;         //lncl/(1-cl)的梯度
double ddlncl;        //lncl/(1-cl)的二阶梯度
double r[1000];//=(rand()/32767.0*2.0-1)*PI/2;
double rr1,rr2,rr,rr3,rr4,rr5,rr6,rr7,rr8,rr9;
//=====================初试化相场和温度场==========================   
         
void IninPhaseAndTemp()   
{   
    int i,j;
        double jifen=0,sum;
    double x1,x2,xx1,xx2,yy1,yy2;
    for(x1=0.001;x1<=0.999;)
         {
         x2=x1+0.001;
     xx1=x1*x1*x1*(10.0-15.0*x1+6.0*x1*x1);
     xx2=x2*x2*x2*(10.0-15.0*x2+6.0*x2*x2);
     yy1=xx1*(1.0-xx1)/(x1*(1.0-x1)*((1.0-xx1)*cle*(1.0-cle)+xx1*cse*(1.0-cse)));
     yy2=xx2*(1.0-xx2)/(x2*(1.0-x2)*((1.0-xx2)*cle*(1.0-cle)+xx2*cse*(1.0-cse)));
         jifen=jifen+(yy1+yy2)*0.001/2;
         x1=x1+0.001;
         }
         sum=jifen*8.314*T*(cle-cse)*(cle-cse)/vm;
         M=jmn*Dl*sqrt(2.0*W)/(sum*eip0*eip0*eip0);
     printf(" %f   %f   %f   %f \n", jmn,  M, Dl, W);

        for(i=0;i<MaxGrid+1;i++)
        {
                for(j=0;j<MaxGrid+1;j++)
                {       
                        for(k=0;k<MaxGrid+1;k++)
                {       

                       
                        if(k<5)//////////if(((i-100)*(i-100)+(j-100)*(j-100)+(k-100)*(k-100))<100)
                        {
                                p[j][k]=1.0;
                            u[j][k]=T;
                                c[j][k]=0.02;
                                newq[j][k]=1.0;
                                newc[j][k]=0.02;
                               
            }
                        else
                        {
                                p[j][k]=0.0;                       
                            u[j][k]=T;
                            c[j][k]=0.02;
                                newq[j][k]=0.0;
                                newc[j][k]=0.02;
                               
                        }
///////////////////////////////////////////////////////////////////////////
                //printf("  %2.3f   \n", p[j][k]);
               
                }
        }
        }
////////////////////////////////////////////////////////////////////
       
///////////////////////////////////////////////////////////////////////
}

//=====================求解相场==========================

void Phase()   
{       
        double px,py,pz,pxx,pxy,pxz,pyy,pyz,pzz,uxx,uyy,angle,angx,angy,angz,angley,anglez,angle1,angle1x,angle1z,angle2,angle2x,angle2y,dajiao,dajiao1,dajiao2;
        double ghyz,ghxz,ghxy,ghxyz,p2x,p2y,p2z;
        double cdj,cdj1,cdj2,sdj,sdj1,sdj2;
           double r;
   
    double gp1,hp1;
    double D1,D2,D3,D4,D5,D6,D7;
        double A,aa,bb;
        double fclcl,fcscs,fccw,fcce,fccs,fccn,fccu,fccd,xj1,xj2,xj3,xj4,xj5,xj6,dc;
        double E12,E11,E13,E1,E2,E3,fq,dq,E12x,E12y,E12z,E1x,E1y,E1z;
        ///////计算相场和温度场///////
   
     
        /*for(i=2;i<MaxGrid+2;i++)
        {
                for(j=2;j<MaxGrid+2;j++)
                {     
                        r[j]=rand()/32767.0*2.0-1;               
                }
        }*/
   
   
  for(i=0;i<MaxGrid+1;i++)
  {
   for(j=0;j<MaxGrid+1;j++)
     {
         for(k=0;k<MaxGrid+1;k++)
     {
   
   
      p[j][k]=newq[j][k];////////***********
          

          c[j][k]=newc[j][k];






     hp=p[j][k]*p[j][k]*p[j][k]*(6*p[j][k]*p[j][k]-15*p[j][k]+10);
       
     
           A=cle*(1.0-cse)/cse/(1.0-cle);
   


         if(p[j][k]<0.001)
         {
         cl[j][k]=c[j][k];
         cs[j][k]=cl[j][k]/(A+(1.0-A)*cl[j][k]);
         }
         else if(p[j][k]>0.999)
         {
                 cs[j][k]=c[j][k];
                 cl[j][k]=A*cs[j][k]/(1.0+(A-1.0)*cs[j][k]);
         }
         else
         {
                 aa=hp*(1.0-A);
                 bb=hp+(1.0-hp)*A+c[j][k]*(1.0-A);
                 cs[j][k]=(bb-sqrt(bb*bb-4.0*aa*c[j][k]))/(2.0*aa);
                 cl[j][k]=(c[j][k]-hp*cs[j][k])/(1.0-hp);
         }
     

     fclcl=8.314*u[j][k]/(cl[j][k]*(1.0-cl[j][k])*vm);
         fcscs=8.314*u[j][k]/(cs[j][k]*(1.0-cs[j][k])*vm);
         fc[j][k]=8.314*u[j][k]*log(cl[j][k]/(1.0-cl[j][k]))/vm;
         fcc[j][k]=fclcl*fcscs/((1.0-hp)*fcscs+hp*fclcl);
   }
  }

  }




  for(i=1;i<MaxGrid;i++)
  {
   for(j=1;j<MaxGrid;j++)
     {
    for(k=1;k<MaxGrid;k++)
     {
     px=(p[i+1][j][k]-p[i-1][j][k])/(2*deltx);

         py=(p[j+1][k]-p[j-1][k])/(2*delty);

         pz=(p[j][k+1]-p[j][k-1])/(2*deltz);
       
         pxy=(p[i+1][j+1][k]-p[i-1][j+1][k]-p[i+1][j-1][k]+p[i-1][j-1][k])/(4*deltx*delty);

         pxz=(p[i+1][j][k+1]-p[i-1][j][k+1]-p[i+1][j][k-1]+p[i-1][j][k-1])/(4*deltx*deltz);

         pyz=(p[j+1][k+1]-p[j+1][k-1]-p[j-1][k+1]+p[j-1][k-1])/(4*deltz*delty);
       
     pxx=(p[i+1][j][k]-2*p[j][k]+p[i-1][j][k])/(deltx*deltx);       
     
     pyy=(p[j+1][k]-2*p[j][k]+p[j-1][k])/(delty*delty);

         pzz=(p[j][k+1]-2*p[j][k]+p[j][k-1])/(deltz*deltz);

     uxx=(u[i+1][j][k]-2*u[j][k]+u[i-1][j][k])/(deltx*deltx);       
       
     uyy=(u[j+1][k]-2*u[j][k]+u[j-1][k])/(delty*delty);
  
     hp1=(30*p[j][k]*p[j][k]*p[j][k]*p[j][k]-60*p[j][k]*p[j][k]*p[j][k]+30*p[j][k]*p[j][k]);
              //hp的一阶导数
      
         gp=p[j][k]*p[j][k]*p[j][k]*p[j][k]-2*p[j][k]*p[j][k]*p[j][k]+p[j][k]*p[j][k];

         gp1=(4*p[j][k]*p[j][k]*p[j][k]-6*p[j][k]*p[j][k]+2*p[j][k]);     //gp的一阶导数

     r=rand()/32767.0*2.0-1;
         
         ghyz=sqrt(py*py+pz*pz);
         ghxz=sqrt(px*px+pz*pz);
         ghxy=sqrt(px*px+py*py);
         ghxyz=sqrt(px*px+py*py+pz*pz);

         p2x=px*px;
         p2y=py*py;
         p2z=pz*pz;

         //printf("%d  %d  %d  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f \n",i,j,k, py,px,pz,pxx,pyy,pzz,pxy,pxz,pyz,ghyz);
     //printf("%f \n", r);
         
/////////////////////////////////////////////////////////////////////////////////                 
     

         
               

     
        if(px!=0)
        {
                angle=atan(ghyz/px);                                         
        }
    else if(ghyz==0)
        {
          angle=PI/4;
        }
        else
        {
          angle=PI/2;
        }         

        if(ghyz==0)
   {
           if(px==0)
         {
         angx=1;
         angy=1;
         angz=1;
         }
           else
           {
           angx=0;    //////pxy/px;
           angy=0;
           angz=0;
          
           }
   }
   else
   {
         angx=(px*py*pxy+px*pz*pxz-pxx*(py*py+pz*pz))/((px*px+py*py+pz*pz)*sqrt(py*py+pz*pz));
     angy=(px*py*pyy+px*pz*pyz-pxy*(py*py+pz*pz))/((px*px+py*py+pz*pz)*sqrt(py*py+pz*pz));
         angz=(px*py*pyz+px*pz*pzz-pxz*(py*py+pz*pz))/((px*px+py*py+pz*pz)*sqrt(py*py+pz*pz));
   }

        if(py!=0)
        {
                angle1=atan(ghxz/py);                                         
        }
    else if(ghxz==0)
        {
          angle1=PI/4;
        }
        else
        {
          angle1=PI/2;
        }         
        if(ghxz==0)
   {
       
           if(py==0)
         {
         angle1x=1;
         angley=1;
         angle1z=1;
         }
           else
           {
           angle1x=0;
           angley=0;          //////pxy/py;;
           angle1z=0;
          
           }
         
   }
   else
   {
         angle1x=(px*py*pxx+py*pz*pxz-pxy*(px*px+pz*pz))/((px*px+py*py+pz*pz)*sqrt(px*px+pz*pz));
     angley=(px*py*pxy+py*pz*pyz-pyy*(px*px+pz*pz))/((px*px+py*py+pz*pz)*sqrt(px*px+pz*pz));
         angle1z=(px*py*pxz+py*pz*pzz-pyz*(px*px+pz*pz))/((px*px+py*py+pz*pz)*sqrt(px*px+pz*pz));
   }
        if(pz!=0)
        {
                angle2=atan(ghxy/pz);                                         
        }
    else if(ghxy==0)
        {
          angle2=PI/4;
        }
        else
        {
          angle2=PI/2;
        }         
        if(ghxy==0)
   {
       
          if(pz==0)
         {
         angle2x=1;
         angle2y=1;
         anglez=1;
         }
           else
           {
           angle2x=0;
           angle2y=0;
           anglez=0;     //////pxz/pz;
          
           }
         
   }
   else
   {
         angle2x=(px*pz*pxx+py*pz*pxy-pxz*(px*px+py*py))/((px*px+py*py+pz*pz)*sqrt(px*px+py*py));
         angle2y=(px*pz*pxy+py*pz*pyy-pyz*(px*px+py*py))/((px*px+py*py+pz*pz)*sqrt(px*px+py*py));
     anglez=(px*pz*pxz+py*pz*pyz-pzz*(px*px+py*py))/((px*px+py*py+pz*pz)*sqrt(px*px+py*py));
         
   }

         
          dajiao=ka*angle;
          cdj=cos(dajiao);
          sdj=sin(dajiao);
          dajiao1=ka*angle1;
          cdj1=cos(dajiao1);
          sdj1=sin(dajiao1);
          dajiao2=ka*angle2;
          cdj2=cos(dajiao2);
          sdj2=sin(dajiao2);
/////////////////////////////////////////////////////////////////////////////////
         if(p[j][k]<0.999)
         {
      
          eipx=eip0*(1+v*cdj);
          eipx1=-eip0*ka*v*sdj;
          eipx2=-eip0*ka*ka*v*cdj;

          eipy=eip0*(1+v*cdj1);
          eipy1=-eip0*ka*v*sdj1;
          eipy2=-eip0*ka*ka*v*cdj1;
         
          eipz=eip0*(1+v*cdj2);
          eipz1=-eip0*ka*v*sdj2;
          eipz2=-eip0*ka*ka*v*cdj2;

          eip=(eipx+eipy+eipz)/3;
          eip1=(eipx1+eipy1+eipz1)/3;
          eip2=(eipx2+eipy2+eipz2)/3;

          E11=(eipx*eipx+eipy*eipy+eipz*eipz)*(pxx+pyy+pzz)/3;

      //E11=eip*eip*(pxx+pyy+pzz);

          E12x=2*eipx*eipx1*(px*angx+py*angy+pz*angz);
          E12y=2*eipy*eipy1*(px*angle1x+py*angley+pz*angle1z);
          E12z=2*eipz*eipz1*(px*angle2x+py*angle2y+pz*anglez);
          E12=(E12x+E12y+E12z)/3;

          E13=angx*ghyz+angley*ghxz+anglez*ghxy;

          E1x=(eipx1*eipx1+eipx*eipx2);
          E1y=(eipy1*eipy1+eipy*eipy2);
          E1z=(eipz1*eipz1+eip*eipz2);

          E1=(E1x+E1y+E1z)/3;
               
          //E1=(eip1*eip1+eip*eip2);
                 
                 
               


         fq=-8.314*u[j][k]*hp1*log((1.0-cse)*(1.0-cl[j][k])/((1.0-cle)*(1.0-cs[j][k])))/vm+W*gp1;
         dq=delt*M*(E11+E12-E1*E13-fq);////我把E2改成-号了

         //printf("%d  %d   %d    %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f  %2.3f \n",i,j,k, angx,angy,angz,angley,anglez,E11,E12,E13,E1,fq);
     /////////////////////////////////////////////////////////  
          if(p[j][k]>0.13 && p[j][k]<0.15)
          {dq=dq+16*gp*r*0.8;}////////**********
          else
          {dq=dq;}
     /////////////////////////////////////////



     newq[j][k]=p[j][k]+dq;
         }
         else
        {
                 newq[j][k]=p[j][k];
         }
         
     if(newq[j][k]<0.0)
          {
                  newq[j][k]=0.0;
          }


       //////////////溶质扩散方程//////////////////////////////


          D1=Ds;
                  D2=Ds;
                  D3=Ds;  
          D4=Ds;
          D5=Ds;
                  D6=Ds;
          D7=Ds;
             



      if(p[j][k]<0.9)

      {
                  D1=Dl;
          }
      
      if(p[i-1][j][k]<0.9)
          {
                  D2=Dl;
          }
          if(p[i+1][j][k]<0.9)
          {
                  D3=Dl;
          }
          if(p[j+1][k]<0.9)
          {
                  D4=Dl;
          }
          if(p[j-1][k]<0.9)
          {
                  D5=Dl;
          }
           if(p[j][k+1]<0.9)
          {
                  D6=Dl;
          }
          if(p[j][k-1]<0.9)
          {
                  D7=Dl;
          }
          
          fccw=2.0/(fcc[i-1][j][k]/D2+fcc[j][k]/D1);
          fcce=2.0/(fcc[i+1][j][k]/D3+fcc[j][k]/D1);
          fccs=2.0/(fcc[j+1][k]/D4+fcc[j][k]/D1);
          fccn=2.0/(fcc[j-1][k]/D5+fcc[j][k]/D1);
          fccu=2.0/(fcc[j][k+1]/D6+fcc[j][k]/D1);
          fccd=2.0/(fcc[j][k-1]/D7+fcc[j][k]/D1);

          xj1=(fc[i-1][j][k]-fc[j][k])*fccw/deltx;
          xj2=(fc[i+1][j][k]-fc[j][k])*fcce/deltx;
          xj3=(fc[j+1][k]-fc[j][k])*fccs/delty;
          xj4=(fc[j-1][k]-fc[j][k])*fccn/delty;
          xj5=(fc[j][k+1]-fc[j][k])*fccu/deltz;
          xj6=(fc[j][k-1]-fc[j][k])*fccd/deltz;
          dc=(xj1+xj2)/deltx+(xj3+xj4)/delty+(xj5+xj6)/deltz;
          newc[j][k]=c[j][k]+dc*delt;
         
          

///printf("%d  %d  %f  %f   %f \n", i,j,p[j],c[j],M);

////////////////////////////////////////////////////////////////////////////////////
         
                 
///////////////////////////////////////////////////////////////////////////////////////////////          
        }
          

  }
}  


  ///////////////设置边界条件//////////
        for(i=0;i<MaxGrid+1;i++)
        {
                for(j=0;j<MaxGrid+1;j++)
        {
                newq[0][j]=newq[1][j];
            
        newq[MaxGrid][j]=newq[MaxGrid-1][j];
        
         
                newq[0][j]=newq[1][j];
            
        newq[MaxGrid][j]=newq[MaxGrid-1][j];
        
               
                newq[j][0]=newq[j][1];
          
        newq[j][MaxGrid]=newq[j][MaxGrid-1];
        
               
        

        u[0][j]=u[1][j];
          
        u[MaxGrid][j]=u[MaxGrid-1][j];
         
                u[0][j]=u[1][j];
            
        u[MaxGrid][j]=u[MaxGrid-1][j];
                      
                u[j][0]=u[j][1];
          
        u[j][MaxGrid]=u[j][MaxGrid-1];
      

               

        newc[0][j]=newc[1][j];
            
        newc[MaxGrid][j]=newc[MaxGrid-1][j];
         
                newc[0][j]=newc[1][j];
            
        newc[MaxGrid][j]=newc[MaxGrid-1][j];
       
                newc[j][0]=newc[j][1];
            
        newc[j][MaxGrid]=newc[j][MaxGrid-1];
      
               
               
               

        f[0][j]=f[1][j];
            
        f[MaxGrid][j]=f[MaxGrid-1][j];
         
                f[0][j]=f[1][j];
            
        f[MaxGrid][j]=f[MaxGrid-1][j];
        
               
                f[j][0]=f[j][1];
            
        f[j][MaxGrid]=f[j][MaxGrid-1];

        
      
                }  

        }
}
//////////////////输出温度场和相场///////////////////

void output_tempresult()
{
        int i,j,k;  
        char filename[30];   
        FILE *fp;

    if(num%10==0)
        {
        sprintf(filename,"%d.plt",num);
    fp=fopen(filename,"w");
   
        fprintf(fp,"%s","title='pu.txt'\n");
        fprintf(fp,"%s","variables= i , j , k ,  p ,   c\n");
        fprintf(fp,"%s","zone I=200,J=200,K=200,F=point\n");  
       
       
                for(i=1;i<MaxGrid+1;i++)
                {
                        for(j=1;j<MaxGrid+1;j++)
                        {
                                for(k=1;k<MaxGrid+1;k++)
                        {
                          fprintf(fp,"%d   %d   %d  ",i,j,k );
                          fprintf(fp,"  %3.2f      %3.4f\n",newq[j][k],newc[j][k]);
              
                        }
                       
                       
                   }

                  
                }
            fprintf(fp,"\n");       
            fclose(fp);
        }
}


int main(int argc, char* argv[])
{
        t=0;
        num=0;
       
//////////////////////////////////////////
printf("this programe is running\n");
//for(i=2;i<MaxGrid+2;i++)
//{
// for(j=2;j<MaxGrid+2;j++)
//        {
     
//     printf("MaxGrid=%1.3d MaxGrid[j]=%1.3d p[j]=%1.3f u[j]=%1.3f c[j]=%1.3f\n",i,j, p[j],u[j],c[j]);
            IninPhaseAndTemp();
//        }
//}
/////////////////////////////////////////   
   while(num<20001)
   {
         Phase();  

    output_tempresult();

    t=t+delt;  
     
        //printf("%d  %1.2f %1.2f %1.2f %1.2f %1.2f %1.2f %1.2f %1.2f %1.2f %1.2f \n", num,rr,rr1,rr2,rr3,rr4,rr5,rr6,rr7,rr8,rr9);
    printf("%d   \n", num);



        num++;
       
   }  
    return 0;
}

MMMIX

回复 3# 长刺的耗子


    编辑下帖子，把表情禁止掉