时间戳

首页 > 解决方案 > 分配数组值时出现分段错误

c - 分配数组值时出现分段错误

问题描述

void particle_initialize()
{
  int i;
    for (i = 0; i<particlenumber; i++) 
{
Pu[i]=0;     
} 
     return; } 

void forceon_particle()
{
for (i = 0; i<particlenumber; i++) 
{
FDx[i]=(U[i]-Pu[i])/taup;
FDy[i]=(V[i]-Pv[i])/taup;
}
return; }

void particle_motion()
{
int i;
for (i = 0; i<particlenumber; i++) 
{
Pu[i]=Pu[i]+FDx[i]*tp;
Px[i]=Px[i]+Pu[i]*tp;}
return; 
}

大家好，我无法理解为什么我在函数particle_motion 中出现分割错误。

只有当我尝试分配时才会出现错误Px[i]=Px[i]+Pu[i]*tp; Py[i]=Py[i]+Pv[i]*tp;。请帮忙。

/*Developing flow in a Two Dimensional channel*/

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#define     NX      101
#define     NY      101
#define     Q       9
#define     T       300001
#define     u_in            0.01
#define     omega            1
#define        kvisc            1/6
#define     W0      4.0/9.0
#define     W1      1.0/9.0
#define     W2      1.0/36.0
#define     rho0        1.0
#define particlenumber 10000
#define dp                   0.01
#define rhop                 2
#define tp                    1/omega                            
#define pi 3.14


int         cx          [Q]={0,1,0,-1,0,1,-1,-1,1};
int         cy          [Q]={0,0,1,0,-1,1,1,-1,-1};
double      t           [Q]={W0,W1,W1,W1,W1,W2,W2,W2,W2};
int         next_x      [NX][Q];
int         next_y      [NY][Q];

double      f1      [NX][NY][Q];
double      f2      [NX][NY][Q];
double      ux      [NX][NY];
double      uy      [NX][NY];
double      uxeq        [NX][NY];//eq velocity for effect of particle
double      uyeq        [NX][NY];
double      ustar       [NX][NY];//forcing term added to ueq
double      vstar       [NX][NY];
double      rho     [NX][NY];
double          Px              [particlenumber];//location of particles
double          Py              [particlenumber];
double          Pu              [particlenumber];//velocity of particles
double          Pv              [particlenumber];

double          FDx              [particlenumber];// drag force
double          FDy              [particlenumber];
double          U              [particlenumber];//velocity of particle
double          V              [particlenumber];
double          FPx             [particlenumber];//force on fluid due to particle
double          FPy             [particlenumber];
double          uforcing             [particlenumber];
double          rhoatpar           [particlenumber];
double          vforcing             [particlenumber];
void        tables();
void        initialise();
void        cal_obs();
void        streaming();
void        collision();
void        density_data();
void particle_initialize();
void forceon_particle();
void intrapolation_velocity();
void particle_motion();
void forceon_fluid();
void intrapolation_density();
void reverse_intrapolation();


int main()
{
    int count=0, k;

    tables();

    initialise();

particle_initialize();

    for(k=0;k<T;k++)
    {
        //printf("k=%d\t",k);
        cal_obs();
        collision();
        streaming();
        if (k % 10000 ==0)
        {
            density_data(count, k);
        }

intrapolation_velocity();
forceon_particle();
intrapolation_density();
forceon_fluid();
reverse_intrapolation();
particle_motion();

    }

    return 0;
}


void tables()
{
    int x,y,i;

    for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
            for (i=0; i<Q; i++)
            {
                next_x[x][i] = x + cx[i];

                    if (next_x[x][i]>NX-1)
                        next_x[x][i] = 0;
                    if (next_x[x][i]<0)
                        next_x[x][i] = NX-1;

                next_y[y][i] = y + cy[i];

                    if (next_y[y][i]>NY-1)
                        next_y[y][i] = 0;
                    if (next_y[y][i]<0)
                        next_y[y][i] = NY-1;
            }
        }
    }

    return;
}

void initialise()
{
    int x,y,i;

    for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
ustar[x][y]=0;//initialize forcing term zero
vstar[x][y]=0;
            for (i=0; i<Q; i++)
            {

                f1[x][y][i]=rho0*t[i];

            }
        }
    }

    return;
}

void cal_obs()
{
    int x,y,i;
    double dense;

    for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
            dense=ux[x][y]=uy[x][y]=0.0;

            for (i=0; i<Q; i++)
            {
                dense     += f1[x][y][i];
                ux[x][y]  += f1[x][y][i]*cx[i];
                uy[x][y]  += f1[x][y][i]*cy[i];
            }

            if(dense != 0.0)
            {
                rho[x][y] = dense;
                ux[x][y]  = ux[x][y]/rho[x][y];
                uy[x][y]  = uy[x][y]/rho[x][y];
            }
            else
            {
                rho[x][y] = 0.0;
                ux[x][y]  = 0.0;
                uy[x][y]  = 0.0;
            }
        }
    }

     for(x=0; x<NX; x++)
    {
        ux[x][0]=u_in;
ux[x][NY-1]=-u_in;
}
    return;
}


void collision()
{
    int x,y,i;
    double udotc,u2,feq;

    for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
            for (i=0; i<Q; i++)
            {
uxeq[x][y]=ustar[x][y]+ux[x][y];
uyeq[x][y]=vstar[x][y]+uy[x][y];
                udotc = uxeq[x][y]*cx[i]+uyeq[x][y]*cy[i];
                u2=uxeq[x][y]*uxeq[x][y] + uyeq[x][y]*uyeq[x][y];

                feq = t[i]*rho[x][y]*(1+3*udotc+4.5*udotc*udotc-1.5*u2);

                f1[x][y][i]=f1[x][y][i]*(1-omega)+feq*omega;
            }
        }
    }

    return;
}


void streaming()
{
    int x,y,i,newx,newy;
    double rho_w;

    for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
            for (i=0; i<Q; i++)
            {
                newx=next_x [x][i];
                newy=next_y[y][i];
                f2[newx][newy][i] = f1[x][y][i];
            }
        }
    }


    for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
            for (i=0; i<Q; i++)
            {
                f1[x][y][i]=f2[x][y][i];
            }
        }
    }

return;
}

void density_data(int data_counter, int id)
{
     int     x, y;

    FILE    *f3;
   char    phase   [100];


/*  Data for Tecplot*/
 sprintf ( phase,"phaseVelocity2D_%d.csv",id);
    f3 = fopen(phase, "w");
   fprintf(f3," X,Y,Rho,U,V\n");
    for (y =0; y < NY; y++)
    {   for (x = 0; x < NX; x++)
            fprintf(f3,"%d%c%d%c%lf%c%lf%c%lf\n",x,',',y,',',rho[x][y],',',ux[x][y],',',uy[x][y]);
        
    }
   
    fflush(f3);
    fclose(f3);
int i;
    FILE    *f4;
   char    particle   [100];

 sprintf ( particle,"particle_%d.csv",id);
    f4 = fopen(particle, "w");
   fprintf(f4," Id,X,Y\n");
    for (i =0; i < particlenumber; i++)
    {   
            fprintf(f4,"%d%c%lf%c%lf\n",i,',',Px[i],',',Py[i]);

    }
   
    fflush(f4);
    fclose(f4);
      return;
}

void particle_initialize()
{
  int i;
srand(time(0));
float randomx;
    for (i = 0; i<particlenumber; i++) 
{
randomx=(float)rand()/RAND_MAX;
Px[i]=randomx*NX;
Pu[i]=0;        
}
   

float randomy;
    for (i = 0; i<particlenumber; i++) 
{
randomy=(float)rand()/RAND_MAX;
Py[i]=randomy*NY;
Pv[i]=0;
        
}
    return; 
} 


void forceon_particle()
{
int i;
double taup;
taup=(rhop*dp*dp)/(18*kvisc*rho0);
for (i = 0; i<particlenumber; i++) 
{
FDx[i]=(U[i]-Pu[i])/taup;
FDy[i]=(V[i]-Pv[i])/taup;
}
return; 

}

void intrapolation_velocity()
{
int i;
int node1x, node1y, node2x, node2y,node3x, node3y,node4x, node4y;
double w1,w2,w3,w4;

for (i = 0; i<particlenumber; i++) 
{

node1x=trunc(Px[i]);
node1y=trunc(Py[i]);
node2x=node1x;
node2y=node1y+1;
node3x=node1x+1;
node3y=node1y+1;
node4x=node1x+1;
node4y=node1y;
w1=(Px[i]-node1x)*(Px[i]-node1x)+(Py[i]-node1y)*(Py[i]-node1y);
w1=sqrt(w1);
w2=(Px[i]-node2x)*(Px[i]-node2x)+(Py[i]-node2y)*(Py[i]-node2y);
w2=sqrt(w2);
w3=(Px[i]-node3x)*(Px[i]-node3x)+(Py[i]-node3y)*(Py[i]-node3y);
w3=sqrt(w3);
w4=(Px[i]-node4x)*(Px[i]-node4x)+(Py[i]-node4y)*(Py[i]-node4y);
w4=sqrt(w4);
U[i]=w1*ux[node1x][node1y]+w2*ux[node2x][node2y]+w3*ux[node3x][node3y]+w4*ux[node4x][node4y];
U[i]=U[i]/(w1+w2+w3+w4);
V[i]=w1*uy[node1x][node1y]+w2*uy[node2x][node2y]+w3*uy[node3x][node3y]+w4*uy[node4x][node4y];
V[i]=V[i]/(w1+w2+w3+w4);

}
return; 
}

void particle_motion()
{

int i;
double A,B;

for (i = 0; i<particlenumber; i++) 
{
Pu[i]=Pu[i]+FDx[i]*tp;//Pu[i] added so that it adds previous velocity
Pv[i]=Pv[i]+FDy[i]*tp;
}

for (i = 0; i<particlenumber; i++) 
{

Px[i]=Px[i]+Pu[i]*tp;
Py[i]=Py[i]+Pv[i]*tp;


//if particle goes out of bound we r introducing periodic BC

if(Px[i]<0||Px[i]>(NX-1))
{
if(Px[i]<0)
{Px[i]=Px[i]+NX-1;}
else if(Px[i]>(NX-1))
{Px[i]=Px[i]-NX+1;}
}


if(Py[i]<0||Py[i]>(NY-1))
{
if(Py[i]<0)
{Py[i]=Py[i]+NY-1;}
else if(Py[i]>(NY-1))
{Py[i]=Py[i]-NY+1;}
}
}

return; 
}

void forceon_fluid()
{
int i;
double taup;
long double mp;
mp=(1/6)*(pi)*rhop*dp*dp*dp;
taup=(rhop*dp*dp)/(18*kvisc*rho0);
for (i = 0; i<particlenumber; i++) 
{
FPx[i]=mp*(U[i]-Pu[i]);
FPx[i]=FPx[i]/taup;
FPy[i]=mp*(V[i]-Pv[i]);
FPy[i]=FPy[i]/taup;
uforcing[i]=(FPx[i])/(rhoatpar[i]*omega);
vforcing[i]=(FPy[i])/(rhoatpar[i]*omega);
}
return;
}


void intrapolation_density()
{
int i;
int node1x, node1y, node2x, node2y,node3x, node3y,node4x, node4y;
double w1,w2,w3,w4;
for (i = 0; i<particlenumber; i++) 
{
node1x=trunc(Px[i]);
node1y=trunc(Py[i]);
node2x=node1x;
node2y=node1y+1;
node3x=node1x+1;
node3y=node1y+1;
node4x=node1x+1;
node4y=node1y;
w1=(Px[i]-node1x)*(Px[i]-node1x)+(Py[i]-node1y)*(Py[i]-node1y);
w1=sqrt(w1);
w2=(Px[i]-node2x)*(Px[i]-node2x)+(Py[i]-node2y)*(Py[i]-node2y);
w2=sqrt(w2);
w3=(Px[i]-node3x)*(Px[i]-node3x)+(Py[i]-node3y)*(Py[i]-node3y);
w3=sqrt(w3);
w4=(Px[i]-node4x)*(Px[i]-node4x)+(Py[i]-node4y)*(Py[i]-node4y);
w4=sqrt(w4);
rhoatpar[i]=w1*rho[node1x][node1y]+w2*rho[node2x][node2y]+w3*rho[node3x][node3y]+w4*rho[node4x][node4y];
rhoatpar[i]=rhoatpar[i]/(w1+w2+w3+w4);
}
return; 
}

void reverse_intrapolation()
{
int i;
int node1x, node1y, node2x, node2y,node3x, node3y,node4x, node4y;
double w1,w2,w3,w4;
int x,y;

for(x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
ustar[x][y]=0;//initialize forcing term zero
vstar[x][y]=0;
}
}

for (i = 0; i<particlenumber; i++) 
{
node1x=trunc(Px[i]);
node1y=trunc(Py[i]);
node2x=node1x;
node2y=node1y+1;
node3x=node1x+1;
node3y=node1y+1;
node4x=node1x+1;
node4y=node1y;
w1=(Px[i]-node1x)*(Px[i]-node1x)+(Py[i]-node1y)*(Py[i]-node1y);
w1=sqrt(w1);
w2=(Px[i]-node2x)*(Px[i]-node2x)+(Py[i]-node2y)*(Py[i]-node2y);
w2=sqrt(w2);
w3=(Px[i]-node3x)*(Px[i]-node3x)+(Py[i]-node3y)*(Py[i]-node3y);
w3=sqrt(w3);
w4=(Px[i]-node4x)*(Px[i]-node4x)+(Py[i]-node4y)*(Py[i]-node4y);
w4=sqrt(w4);

ustar[node1x][node1y]=ustar[node1x][node1y]+(w1*uforcing[i])/(w1+w2+w3+w4);
vstar[node1x][node1y]=vstar[node1x][node1y]+(w1*vforcing[i])/(w1+w2+w3+w4);

ustar[node2x][node2y]=ustar[node2x][node2y]+(w1*uforcing[i])/(w1+w2+w3+w4);
vstar[node2x][node2y]=vstar[node2x][node2y]+(w1*vforcing[i])/(w1+w2+w3+w4);

ustar[node3x][node3y]=ustar[node3x][node3y]+(w1*uforcing[i])/(w1+w2+w3+w4);
vstar[node3x][node3y]=vstar[node3x][node3y]+(w1*vforcing[i])/(w1+w2+w3+w4);

ustar[node4x][node4y]=ustar[node4x][node4y]+(w1*uforcing[i])/(w1+w2+w3+w4);
vstar[node4x][node4y]=vstar[node4x][node4y]+(w1*vforcing[i])/(w1+w2+w3+w4);

}
return; 
}

我已经添加了整个代码，因为错误可能是由于其他函数也可能访问数组。

标签： c

解决方案

问题在于这个声明。

V[i] = w1* uy[node1x][node1y] + w2* uy[node2x][node2y] + w3 * uy[node3x][node3y] + w4*uy[node4x][node4y];

这里，数组是uy[101][101], ux[101][101]。查看 和 的node1x值node1y：

node1x=trunc(Px[i]);
node1y=trunc(Py[i]);

查看 和 的Px[i]值Py[i]：

Px[i]=randomx*NX; // NX=101. randomx ranges from 0 to 1. randomx=(float)rand()/RAND_MAX
Py[i]=randomy*NY; // NY=101. randomy ranges from 0 to 1. randomy=(float)rand()/RAND_MAX

我想在某个时间点randomx = 1，那么Px[i] = 101，node1x = 101那么你正在访问uy[101][node1y]。这会导致分段错误，因为数组是uy[101][101]，并且您无法访问uy[101][node1y]元素。

推荐阅读