分子流れの数値シミュレーションをDSMC法を用いて行っています。
// maximum collision number method ////////
の
g[k]=c[n][k]-c[l][k];
のところでエラーが出ます。「…でハンドルされていない例外が発生しました:… を読み込み中にアクセス違反が発生しました。」
というやつです。
値を見てみると、
m=(int)(((double)rand()/(double)RAND_MAX)*nm+mcr[j]+1);
でmの値が1201となっており、lcr[m]がlcr[1201]となり(lcrは0~1200の配列で1201は定義されていません)、それがc[n][k],c[l][k]のlやnに代入されることにより発生しているのかと考えております(c[][]も[0~1200][0~3]の配列です)。
ちなみに全粒子数は1200でlcr[]は境界で反射する粒子の数を記憶するための配列なので、そもそも1200になるはずはありません。
mの式がおかしいのでしょうか?乱数がおかしいのでしょうか?型をint型に変換しているのがおかしいのでしょうか?
ちなみに、これは「原子・分子の流れ」(日本機械学会)のfortranで書かれたプログラムをcに自分で書き直したものです。
[2] 環境
[2.1] OS : Windows
[2.2] コンパイラ名 : VC++ 2010Express
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define pi 3.1415927
#define pi2 2.0*pi
#define avog 6.02217*pow(10.0,23.0)
#define boltz 1.38066*pow(10.0,-23.0)
#define wmol 28.013*pow(10.0,-3.0)
#define w wmol/avog
#define cgas boltz/w
#define temp0 273.0
#define u0 1.0
#define hu0 u0/2.0
#define akn 0.05
#define nc 40
#define n0 30
#define nt n0*nc
#define el 1.0
#define hel el/2.0
#define dy el/nc
#define alp 0.2
#define dt0 sqrt(pi)/2.0*alp*akn
#define gmx 5.0/sqrt(2.0)
#define cmx gmx*dt0/(sqrt(8.0)*n0*akn)
int main(void){
int it;
int i;
int j;
int k;
int m;
int n;
int nsum[nc+1];
double sum[4+1][nc+1];
for(j=1; j<=nc; j++){
nsum[j]=0;
for(k=1; k<=4; k++){
sum[k][j]=0.0;
}
}
double y[nt+1];
double y2[nc+1];
double c[nt+1][3+1];
double a;
double b;
i=0;
for(j=1; j<=nc; j++){
y2[j]=-1.0*hel+dy*(j-1);
for(m=1; m<=n0; m++){
i=i+1;
y[i]=y2[j]+dy*((double)rand()/(double)RAND_MAX);
a=sqrt(-1.0*log((double)rand()/(double)RAND_MAX));
b=pi2*((double)rand()/(double)RAND_MAX);
c[i][1]=a*cos(b);
c[i][2]=a*sin(b);
a=sqrt(-1.0*log((double)rand()/(double)RAND_MAX));
b=pi2*((double)rand()/(double)RAND_MAX);
c[i][3]=a*sin(b);
}
}
int itm=6000;
int ismp=5;
int istart=10000+5;
int nsmp=0;
int ip;
double yy[nt+1];
int lcr[nt+1];
double dtres[nt+1];
int iq;
double sgn;
double dt;
double dt1;
int mcr[nt+1];
int ncl[nc+1];
int icell[nt+1];
int nm;
int coll;
int nmax;
int nn;
int l;
double gm;
double g[3+1];
double cs;
double sn;
double rd[3+1];
double cd[3+1];
// advance time ////////
for(it=1; it<=itm; it++){
// free motion ////////
ip=0;
for(i=1; i<=nt; i++){
yy[i]=y[i]+c[i][2]*dt0;
if(fabs(yy[i])<hel) goto skip10;
ip=ip+1;
lcr[ip]=i;
skip10:;
}
// reflection at y=+L/2,-L/2 ////////
for(m=1; m<=ip; m++){
dtres[m]=dt0;
}
skip300:;
iq=0;
for(m=1; m<=ip; m++){
i=lcr[m];
sgn=+1.0;
if(c[i][2]<0.0) sgn=-1.0;
dt1=(sgn*hel-y[i])/c[i][2];
dt=dtres[m]-dt1;
c[i][2]=-1.0*sgn*sqrt(-1.0*log((double)rand()/(double)RAND_MAX));
a=sqrt(-1.0*log((double)rand()/(double)RAND_MAX));
b=pi2*((double)rand()/(double)RAND_MAX);
c[i][1]=a*cos(b)+sgn*hu0;
c[i][3]=a*sin(b);
y[i]=sgn*hel;
yy[i]=y[i]+c[i][2]*dt;
// chek of double, triple, ... reflection ////////
if(fabs(yy[i])<hel) goto skip20;
iq=iq+1;
dtres[iq]=dt;
mcr[iq]=i;
skip20:;
}
if(iq==0) goto skip310;
ip=iq;
for(m=1; m<=iq; m++){
lcr[m]=mcr[m];
}
goto skip300;
skip310:;
for(i=1; i<=nt; i++){
y[i]=yy[i];
}
// indexing ////////
// ncl is number of molecules in a cell ////////
for(j=1; j<=nc; j++){
ncl[j]=0;
mcr[j]=0;
}
for(i=1; i<=nt; i++){
j=(int)((y[i]+hel)/dy+1.0);
// for safety ////////
if(j<1) j=1;
if(j>nc) j=nc;
icell[i]=j;
ncl[j]=ncl[j]+1;
}
n=0;
// mcr[j] is the number of molecules lying before cell j ////////
for(j=1; j<=nc; j++){
mcr[j]=n;
n=n+ncl[j];
ncl[j]=0;
}
for(i=1; i<=nt; i++){
j=icell[i];
ncl[j]=ncl[j]+1;
k=mcr[j]+ncl[j];
lcr[k]=i;
}
// maximum collision number method ////////
for(j=1; j<=nc; j++){
nm=ncl[j];
coll=(int)(cmx*nm*(nm-1));
nmax=coll;
if(((double)rand()/(double)RAND_MAX)<(coll-nmax)) nmax=nmax+1;
for(nn=1; nn<=nmax; nn++){
m=(int)(((double)rand()/(double)RAND_MAX)*nm+mcr[j]+1);
l=lcr[m];
skip320:
m=(int)(((double)rand()/(double)RAND_MAX)*nm+mcr[j]+1);
n=lcr[m];
if(n==1) goto skip320;
gm=0.0;
for(k=1; k<=3; k++){
g[k]=c[n][k]-c[l][k];
gm=gm+g[k]*g[k];
}
gm=sqrt(gm);
a=gm/gmx;
if(((double)rand()/(double)RAND_MAX)>a) goto skip40;
cs=1.0-2.0*((double)rand()/(double)RAND_MAX);
sn=sqrt(1.0-cs*cs);
b=pi2*((double)rand()/(double)RAND_MAX);
rd[1]=sn*cos(b);
rd[2]=sn*sin(b);
rd[3]=cs;
for(k=1; k<=3; k++){
cd[k]=c[l][k]+c[n][k];
}
for(k=1; k<=3; k++){
c[l][k]=0.5*(cd[k]-gm*rd[k]);
c[n][k]=0.5*(cd[k]+gm*rd[k]);
}
skip40:;
}
}
// sampling ////////
if(it<istart) goto skip100;
if(it%ismp!=0) goto skip100;
nsmp=nsmp+1;
for(j=1; j<=nc; j++){
nsum[j]=nsum[j]+ncl[j];
}
for(i=1; i<=nt; i++){
a=0.0;
j=icell[i];
for(k=1; k<=3; k++){
sum[k][j]=sum[k][j]+c[i][k];
a=a+c[i][k]*c[i][k];
}
sum[4][j]=sum[4][j]+a;
}
skip100:;
}
double usq;
double temp;
double density;
double yc;
double u[4+1];
for(j=1; j<=nc; j++){
yc=-1.0*hel+dy*(j-0.5);
for(k=1; k<=4; k++){
u[k]=sum[k][j]/nsum[j];
}
usq=u[1]*u[1]+u[2]*u[2]+u[3]*u[3];
temp=2.0/3.0*(u[4]-usq);
density=nsum[j]/(n0*nsmp);
printf("%f\n%f\n%f\n", usq, temp, density);
}
return 0;
}