DSMC法での粒子衝突におけるエラー

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
grass

DSMC法での粒子衝突におけるエラー

#1

投稿記事 by grass » 5年前

C言語数か月の初心者です。
分子流れの数値シミュレーションを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;
}

アバター
みけCAT
記事: 6246
登録日時: 9年前
住所: 千葉県
連絡を取る:

Re: DSMC法での粒子衝突におけるエラー

#2

投稿記事 by みけCAT » 5年前

原因になっているかはわかりませんが、とりあえず
grass さんが書きました:

コード:

#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
ここの記述において、
例えばw=28.013*pow(10.0,-3.0)/6.02217*pow(10.0,23.0)
cgas=1.38066*pow(10.0,-23.0)/28.013*pow(10.0,-3.0)/6.02217*pow(10.0,23.0)
となりますが、これは意図した動作でしょうか?
w=(28.013*pow(10.0,-3.0))/(6.02217*pow(10.0,23.0))
cgas=(1.38066*pow(10.0,-23.0))/((28.013*pow(10.0,-3.0))/(6.02217*pow(10.0,23.0)))
という値にしたいわけではないでしょうか?
この下の部分のマクロにも、同様の疑問があります。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

grass

Re: DSMC法での粒子衝突におけるエラー

#3

投稿記事 by grass » 5年前

ありがとうございます!

全てのマクロを見直して修正しました。

が、依然同じ箇所で同じエラーがでます…。

grass

Re: DSMC法での粒子衝突におけるエラー

#4

投稿記事 by grass » 5年前

変数の値を調べてみました。

コード:

// 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=((double)rand()/(double)RAND_MAX)*nm+mcr[j]+1;
			l=lcr[m];
			skip320:
			m=((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];
			}
この部分で、j=40のとき、mが1191、lcr[1191]が919の値であるにもかかわらず、l=lcr[m]、つまりl=919が正常に作動せず、lは-858993460となってしまっています。時間進展では、it=1149まで、また上のループでもj=40まで進んでいるので、正常に作動する場合もあるようです。

アバター
へにっくす
記事: 628
登録日時: 7年前
住所: 東京都

Re: DSMC法での粒子衝突におけるエラー

#5

投稿記事 by へにっくす » 5年前

grass さんが書きました:ありがとうございます!

全てのマクロを見直して修正しました。

が、依然同じ箇所で同じエラーがでます…。
たとえ同じエラーでも、
書きなおしたのであれば、再度コードを掲示してください。
でないと本当に直したのか分かりませんよ

また、C言語では単に宣言しただけでは、どういう値が入っているか分かりません。
このことは、

コード:

int a;
printf("a=%d", a);
と値を出してみるだけでもすぐ分かるはずです。
(偶然0となる場合もあるかもしれませんが、明示的に0とは書いていないので駄目です)

コード:

int a = 0;
int b[1000] = {0};
というように必ず初期化しておくことをお勧めします。
オフトピック
もっと言うと、goto文使ってるのが、流れを追いにくくしていますね。
written by へにっくす

grass

Re: DSMC法での粒子衝突におけるエラー

#6

投稿記事 by grass » 5年前

きちんと書いてなくてすみません。
以下のように修正しました。

依然同じ箇所で同じエラーがでています。

コード:

#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=0;
int i=0;
int j=0;
int k=0;
int l=0;
int m=0;
int n=0;
double r=0.0;
int nsum[nc+1]={0};
double sum[4+1][nc+1]={{0.0}};


for(j=1; j<=nc; j++){
	nsum[j]=0;
	for(k=1; k<=4; k++){
		sum[k][j]=0.0;
	}
}


double y[nt+1]={0.0};
double y2[nc+1]={0.0};
double c[nt+1][3+1]={{0.0}};
double a=0.0;
double b=0.0;

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;
		r=(double)rand()/(double)RAND_MAX;
		y[i]=y2[j]+dy*r;
		r=(double)rand()/(double)RAND_MAX;
		a=sqrt(-1.0*log(r));
		r=(double)rand()/(double)RAND_MAX;
		b=pi2*r;
		c[i][1]=a*cos(b);
		c[i][2]=a*sin(b);
		r=(double)rand()/(double)RAND_MAX;
		a=sqrt(-1.0*log(r));
		r=(double)rand()/(double)RAND_MAX;
		b=pi2*r;
		c[i][3]=a*sin(b);
	}
}

int itm=6000;
int ismp=5;
int istart=10000+5;
int nsmp=0;
int ip=0;
double yy[nt+1]={0.0};
int lcr[nt+1]={0};
double dtres[nt+1]={0.0};
int iq=0;
double sgn=0.0;
double dt=0.0;
double dt1=0.0;
int mcr[nt+1]={0};
int ncl[nc+1]={0};
int icell[nt+1]={0};
int nm=0;
int coll=0;
int nmax=0;
int nn=0;
double gm=0.0;
double g[3+1]={0.0};
double cs=0.0;
double sn=0.0;
double rd[3+1]={0.0};
double cd[3+1]={0.0};

// 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;
		r=(double)rand()/(double)RAND_MAX;
		c[i][2]=-1.0*sgn*sqrt(-1.0*log(r));
		r=(double)rand()/(double)RAND_MAX;
		a=sqrt(-1.0*log(r));
		r=(double)rand()/(double)RAND_MAX;
		b=pi2*r;
		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;
		r=(double)rand()/(double)RAND_MAX;
		if(r<(coll-nmax))	nmax=nmax+1;
		
		for(nn=1; nn<=nmax; nn++){
			r=(double)rand()/(double)RAND_MAX;
			m=(int)(r*nm+mcr[j]+1);
			l=lcr[m];
			skip320:
			r=(double)rand()/(double)RAND_MAX;
			m=(int)(r*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;
			r=(double)rand()/(double)RAND_MAX;
			if(r>a)	goto skip40;
			r=(double)rand()/(double)RAND_MAX;
			cs=1.0-2.0*r;
			sn=sqrt(1.0-cs*cs);
			r=(double)rand()/(double)RAND_MAX;
			b=pi2*r;
			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;
}

can110

Re: DSMC法での粒子衝突におけるエラー

#7

投稿記事 by can110 » 5年前

(double)rand()/(double)RAND_MAX が1の場合に範囲外の配列位置を算出しているようですね。
FORTRANとは乱数処理の挙動が異なるのかもしれません。
(double)rand()/(RAND_MAX+1) に変更することで0以上1未満の値が得られ、236行目まではエラー発生しなくなりました。
しかしサンプリング?が行われないため、表示処理にて0除算してしまいます。
219行目の判定は明らかにおかしいと思われるのでistart = 0と修正することにより、以下の結果が得られました。

0.192947
1.021943
1.000000
0.166684
1.029781
1.000000
0.143056
1.030630
1.000000

いかがでしょうか?

grass

Re: DSMC法での粒子衝突におけるエラー

#8

投稿記事 by grass » 5年前

乱数についてご指摘ありがとうございます。
指示の通り、乱数生成をすべて
r=(double)rand()/(RAND_MAX+1);

のように書き直したところ、エラーは発生しなくなりました。
が、計算結果が出力されません。

以下、乱数生成以外にミスを見つけて修正した箇所です。
行数は以下に貼りなおしたコードでの行数です。

59行目:
y2[]と配列となっていましたが、変数でした。→yj=-1.0*hel+dy*(j-1);
78行目:
int itm=6000;→int itm=60000;

この修正を行ったところ、エラーはないのですが、計算がいつまでたっても終わりません。
警告として、警告:式の整数がオーバーフローしました
と出ます。

ご指摘のサンプリングについてですが、
234行目
if(it<istart) goto skip100;
は、ループ回数が10000+5未満の場合、サンプリングを行わずskip100まで飛ばし、
235行目
if(it%ismp!=0) goto skip100;
は、ループ回数がismp=5の倍数でないときもサンプリングを行わずskip100まで飛ばすという意味で作成しております。

istartはサンプリングを開始するループ番号を意味していて、10000+5回目から5回おきに60000回目までサンプリングを行うという意図でコードを作成しました。

一歩前進しましたが、まだうまくいきません。
アドバイスよろしくお願いいたします。

コード:

#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=0;
int i=0;
int j=0;
int k=0;
int l=0;
int m=0;
int n=0;
double r=0.0;
int nsum[nc+1]={0};
double sum[4+1][nc+1]={{0.0}};


for(j=1; j<=nc; j++){
	nsum[j]=0;
	for(k=1; k<=4; k++){
		sum[k][j]=0.0;
	}
}


double y[nt+1]={0.0};
double yj=0.0;
double c[nt+1][3+1]={{0.0}};
double a=0.0;
double b=0.0;

i=0;
for(j=1; j<=nc; j++){
	yj=-1.0*hel+dy*(j-1);
	for(m=1; m<=n0; m++){
		i=i+1;
		r=(double)rand()/(RAND_MAX+1);
		y[i]=yj+dy*r;
		r=(double)rand()/(RAND_MAX+1);
		a=sqrt(-1.0*log(r));
		r=(double)rand()/(RAND_MAX+1);
		b=pi2*r;
		c[i][1]=a*cos(b);
		c[i][2]=a*sin(b);
		r=(double)rand()/(RAND_MAX+1);
		a=sqrt(-1.0*log(r));
		r=(double)rand()/(RAND_MAX+1);
		b=pi2*r;
		c[i][3]=a*sin(b);
	}
}

int itm=60000;
int ismp=5;
int istart=10000+5;
int nsmp=0;
int ip=0;
double yy[nt+1]={0.0};
int lcr[nt+1]={0};
double dtres[nt+1]={0.0};
int iq=0;
double sgn=0.0;
double dt=0.0;
double dt1=0.0;
int mcr[nt+1]={0};
int ncl[nc+1]={0};
int icell[nt+1]={0};
int nm=0;
int coll=0;
int nmax=0;
int nn=0;
double gm=0.0;
double g[3+1]={0.0};
double cs=0.0;
double sn=0.0;
double rd[3+1]={0.0};
double cd[3+1]={0.0};

// 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;
		r=(double)rand()/(RAND_MAX+1);
		c[i][2]=-1.0*sgn*sqrt(-1.0*log(r));
		r=(double)rand()/(RAND_MAX+1);
		a=sqrt(-1.0*log(r));
		r=(double)rand()/(RAND_MAX+1);
		b=pi2*r;
		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;
		r=(double)rand()/(RAND_MAX+1);
		if(r<(coll-nmax))	nmax=nmax+1;
		
		for(nn=1; nn<=nmax; nn++){
			r=(double)rand()/(RAND_MAX+1);
			m=(int)(r*nm+mcr[j]+1);
			l=lcr[m];
			skip320:
			r=(double)rand()/(RAND_MAX+1);
			m=(int)(r*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;
			r=(double)rand()/(RAND_MAX+1);
			if(r>a)	goto skip40;
			r=(double)rand()/(RAND_MAX+1);
			cs=1.0-2.0*r;
			sn=sqrt(1.0-cs*cs);
			r=(double)rand()/(RAND_MAX+1);
			b=pi2*r;
			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;
}

can110

Re: DSMC法での粒子衝突におけるエラー

#9

投稿記事 by can110 » 5年前

サンプリング開始条件について了解しました。

まずは100行目付近に以下のような出力文を挿入してどこまで処理が進んでいるのか確認してみましょう。

コード:

printf( "it=%d\n", id);
処理が止まっている場所が特定できると思います。

次に以下のように条件判定を追加しブレークポイントを貼り、デバッグ実行にてステップ実行してみます。

コード:

if( it == ?????){ // 止まったときの値
    printf( "it=%d\n", id); // この行にブレークポイントを設定する
}
各変数の値を細かく確認してみると、「異常な値」が格納されている変数が見つかります。
なぜ「異常な値」が格納されたのでしょう?
格納元の処理を追っていくと、数学上問題のある演算をしている場所が見つかります。

結論からいうと乱数処理の修正により正常に終了するようになり、以下のような出力が得られました。
0.192500
1.027083
1.000000
0.168236
1.037314
1.000000
(以下略)
乱数処理を修正したため、修正方法によっては異なる結果になりえます。
またこの結果が正しいかどうかはわかりません。

can110

Re: DSMC法での粒子衝突におけるエラー

#10

投稿記事 by can110 » 5年前

100行目付近→107行目付近が正しく、printf中の変数名id→itでした。

コード:

printf( "it=%d\n", it);
失礼しました。

アバター
みけCAT
記事: 6246
登録日時: 9年前
住所: 千葉県
連絡を取る:

Re: DSMC法での粒子衝突におけるエラー

#11

投稿記事 by みけCAT » 5年前

can110 さんが書きました:まずは100行目付近に以下のような出力文を挿入してどこまで処理が進んでいるのか確認してみましょう。

コード:

printf( "it=%d\n", id);
処理が止まっている場所が特定できると思います。
今回は強制終了が発生するため、念のため

コード:

printf( "it=%d\n", it);fflush(stdout);
とした方がいいかもしれません。

【訂正】よく読んでいませんでした。必要ないかもしれません。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

閉鎖

“C言語何でも質問掲示板” へ戻る