/*----------------------------------------*/ /* Riemann problem    */ /*---------------------------------------*/ #include #include #include #define n 100 /* 書き出すメッシュの数 */ #define Length 1.0e0 /* 一辺の長さ */ #define rhoL 1.e0 /* 左側の初期密度 */ #define pL 1.e0 /* 左側の初期圧力 */ #define vL 0.e0 /* 左側の初期速度 */ #define cL sqrt(gamma*pL/rhoL) /* 左側の初期音速 */ #define rhoR 0.125e0 /* 右側の初期密度 */ #define pR 0.1e0 /* 右側の初期圧力 */ #define vR 0.e0 /* 右側の初期速度 */ #define cR sqrt(gamma*pR/rhoR) /* 右側の初期音速 */ #define gamma (4.e0/3.e0)/* 比熱比 */ #define eps 1.e-6 /* この値よりずれが小さくなればOK */ #define pi 3.141592 /* π */ #define tcal 0.2 /* 計算時間 */ #define g1 (gamma-1.e0)/(2.e0*gamma) /* 比熱比から計算したもの */ #define g2 (gamma+1.e0)/(2.e0*gamma) /* 比熱比から計算したもの */ #define g3 (2.e0*gamma)/(gamma-1.e0) /* 比熱比から計算したもの */ #define g4 2.e0/(gamma-1.e0) /* 比熱比から計算したもの */ #define g5 2.e0/(gamma+1.e0) /* 比熱比から計算したもの */ #define g6 (gamma-1.e0)/(gamma+1.e0) /* 比熱比から計算したもの */ #define g7 (gamma-1.e0)/2.e0 /* 比熱比から計算したもの */ #define g8 (gamma-1.e0) /* 比熱比から計算したもの */ void starcal(double *pM, double *vM); /* star 領域の圧力、速度の計算 */ void prefunc(int istat, double *F, double *dF, double p); /* star 領域の圧力を求めるための式を与える */ double pstart(void); /* star 領域の圧力の初期値を与える */ void output(double pM, double vM, double s, double *rho, double *p, double *v,int *ii); /* リーマン問題の解を書き出す */ main(void) { int i,num; double delta_L,dt; double FL,dFL,fR,dFR; double pM,vM,dx,x,s,rho,p,v; int ii; FILE *output_file; output_file=fopen("output.dat","w"); /* 出力ファイルオープン*/ starcal(&pM,&vM); dx=Length/(float)n; for(i=0; isTL){ *rho=rhoL*pow((float)(pM/pL),(float)(1.e0/gamma)); *v=vM; *p=pM; *ii=2; }else{ *v=g5*(cL+g7*vL+s); c=g5*(cL+g7*(vL-s)); *rho=rhoL*pow((c/cL),g4); *p=pL*pow((c/cL),g3); *ii=3; } } }else{ pML=pM/pL; sL=vL-cL*sqrt(g2*pML+g1); if(s<=sL){ *rho=rhoL; *v=vL; *p=pL; *ii=4; }else{ *rho=rhoL*(pML+g6)/(pML*g6+1.e0); *v=vM; *p=pM; *ii=5; } } }else{ if(pM>pR){ pMR=pM/pR; sR=vR+cR*sqrt(g2*pMR+g1); if(s>=sR){ *rho=rhoR; *v=vR; *p=pR; *ii=6; }else{ *rho=rhoR*(pMR+g6)/(pMR*g6+1.e0); *v=vM; *p=pM; *ii=7; } }else{ sHR=vR+cR; if(s>=sHR){ *rho=rhoR; *v=vR; *p=pR; *ii=8; }else{ cMR=cR*pow((pM/pR),g1); sTR=vM+cMR; if(s<=sTR){ *rho=rhoR*pow((pM/pR),(1.e0/gamma)); *v=vM; *p=pM; *ii=9; }else{ *v=g5*(-cR+g7*vR+s); c=g5*(cR-g7*(vR-s)); *rho=rhoR*pow((c/cR),g4); *p=pR*pow((c/cR),g3); *ii=10; } } } } }