fork download
copy 
  1. #include<stdio.h>
  2. #include<math.h>
  3.  
  4.  
  5. int main(){
  6.  
  7. 	int i,j,k,l,m,n;
  8. 	double gamma=1.4,dltt=0.1,dltx=0.1;
  9. 	double q[3][100],qnew[3][100];
  10. 	double E1[100],E2[100],E3[100];
  11. 	double count[100],rho[100],u[100],p[100],ene[100],H[100];
  12. 	double rhoave[99],uave[99],cave[99],Have[99];
  13. 	double detA[99],R[3][3],Ri[3][3];
  14. 	double lambda[3],b1,b2,A[3][3]={0.0};
  15.  
  16. 	FILE *fp;
  17. 	fp=fopen("FDS.csv","w");
  18.  
  19. 	//初期条件//
  20. 	for(j=0;j<51;j++){
  21. 		count[j]=j;
  22. 		rho[j] = 1.0;
  23. 		u[j]   = 0.0;
  24. 		p[j]   = 0.8;
  25. 	}
  26. 	for(j=51;j<100;j++){
  27. 		count[j]=j;
  28. 		rho[j] = 0.125;
  29. 		u[j]   = 0.0;
  30. 		p[j]   = 0.1;
  31. 	}
  32. 	for(j=0;j<100;j++){
  33. 		ene[j] = p[j]/(gamma - 1.0) + 0.5*rho[j]*u[j]*u[j];
  34. 		H[j] = (ene[j]+p[j])/rho[j];
  35. 	}
  36. 	for(j=0;j<100;j++){
  37. 		q[0][j] = rho[j];
  38. 		q[1][j] = rho[j]*u[j];
  39. 		q[2][j] = ene[j];
  40. 	}
  41.  
  42.  
  43.  
  44.  
  45. 	//メインループ//
  46. 	for(m=0;m<10;m++){
  47. 		//平均値//
  48. 		for(j=0;j<99;j++){
  49. 			rhoave[j] = sqrt(rho[j]*rho[j+1]);
  50. 		}
  51. 		for(j=0;j<99;j++){
  52. 			uave[j] = (sqrt(rho[j])*u[j] + sqrt(rho[j+1])*u[j+1])/(sqrt(rho[j]) + sqrt(rho[j+1]));
  53. 		}
  54. 		for(j=0;j<99;j++){
  55. 			Have[j] = (sqrt(rho[j])*H[j] + sqrt(rho[j+1])*H[j+1])/(sqrt(rho[j])+sqrt(rho[j+1]));
  56. 		}
  57. 		for(j=0;j<99;j++){
  58. 			cave[j] = sqrt((gamma - 1.0)*(Have[j] - 0.5*uave[j]*uave[j]));
  59. 		}
  60.  
  61. 		//行列//
  62. 		//定義//
  63. 		for(j=0;j<100;j++){
  64.  
  65. 			lambda[0]=fabs(uave[j] - cave[j]);
  66. 			lambda[1]=fabs(uave[j]);
  67. 			lambda[2]=fabs(uave[j] + cave[j]);
  68.  
  69. 			R[0][0]=1.0;
  70. 			R[0][1]=1.0;
  71. 			R[0][2]=1.0;
  72. 			R[1][0]=uave[j] - cave[j];
  73. 			R[1][1]=uave[j];
  74. 			R[1][2]=uave[j] + cave[j];
  75. 			R[2][0]=Have[j] - uave[j]*cave[j];
  76. 			R[2][1]=0.5*uave[j]*uave[j];
  77. 			R[2][2]=Have[j] + uave[j]*cave[j];
  78.  
  79. 			b1 = 0.5*(gamma - 1.0)*uave[j]*uave[j]/(cave[j]*cave[j]);
  80. 			b2 = (gamma - 1.0)/(cave[j]*cave[j]);
  81.  
  82.  
  83. 			Ri[0][0]=  0.5*(b1 + uave[j]/cave[j]);
  84. 			Ri[0][1]= -0.5*(1.0/cave[j] + b2*uave[j]);
  85. 			Ri[0][2]=  0.5*b2;
  86. 			Ri[1][0]=  1.0 - b1;
  87. 			Ri[1][1]=  b2*uave[j];
  88. 			Ri[1][2]= -b2;
  89. 			Ri[2][0]=  0.5*(b1 - uave[j]/cave[j]);
  90. 			Ri[2][1]=  0.5*(1.0/cave[j] - b2*uave[j]);
  91. 			Ri[2][2]=  0.5*b2;
  92. 			//行列計算//
  93. 			for(l=0;l<3;l++){
  94. 				for(m=0;m<3;m++){
  95. 					A[l][m]=0.0;
  96. 					for(n=0;n<3;n++){
  97. 						A[l][m] += R[l][n]*lambda[n]*Ri[n][m];
  98. 					}
  99. 				}
  100. 			}
  101.  
  102. 			E1[j]=0.5*(q[1][j+1]+q[1][j]
  103. 			-A[0][0]*(q[0][j+1]-q[0][j])
  104. 				-A[0][1]*(q[1][j+1]-q[1][j])
  105. 				-A[0][2]*(q[2][j+1]-q[2][j]));
  106.  
  107. 			E2[j]=0.5*((gamma - 1.0)*(q[2][j+1]+q[2][j] - 0.5*(q[1][j+1]*q[1][j+1]/(q[0][j+1]*q[0][j+1]) + q[1][j]*q[1][j]/(q[0][j]*q[0][j]))) + (q[1][j+1]*q[1][j+1]/q[0][j+1]+q[1][j]*q[1][j]/q[0][j])
  108. 				-A[1][0]*(q[0][j+1]-q[0][j])
  109. 				-A[1][1]*(q[1][j+1]-q[1][j])
  110. 				-A[1][2]*(q[2][j+1]-q[2][j]));
  111.  
  112. 			E3[j]=0.5*((q[2][j+1] + (gamma - 1.0)*(q[2][j+1] - 0.5*q[1][j+1]*q[1][j+1]/q[0][j+1]))*(q[1][j+1]/q[0][j+1]) + (q[2][j] + (gamma - 1.0)*(q[2][j] - 0.5*q[1][j]*q[1][j]/q[0][j]))*(q[1][j]/q[0][j])
  113. 				-A[2][0]*(q[0][j+1]-q[0][j])
  114. 				-A[2][1]*(q[1][j+1]-q[1][j])
  115. 				-A[2][2]*(q[2][j+1]-q[2][j]));
  116.  
  117. 		}
  118. 		//次ステップのQの計算//
  119. 		for(j=1;j<99;j++){
  120. 			qnew[0][j]=q[0][j]-(E1[j]-E1[j-1]);
  121. 			qnew[1][j]=q[1][j]-(E2[j]-E2[j-1]);
  122. 			qnew[2][j]=q[2][j]-(E3[j]-E3[j-1]);
  123. 		}
  124.  
  125. 		//端点//
  126. 		for(i=0;i<3;i++){
  127. 			qnew[i][0]=q[i][0];
  128. 			qnew[i][99]=q[i][99];
  129. 		}
  130.  
  131. 		//qnewの中身をrho,u,pなどに変換//
  132. 		for(j=0;j<100;j++){
  133. 			rho[j]=qnew[0][j];
  134. 			u[j]  =qnew[1][j]/qnew[0][j];
  135. 			p[j]  =(gamma - 1.0)*(qnew[2][j] - 0.5*qnew[1][j]*qnew[1][j]/qnew[0][j]);
  136. 			H[j]  =(qnew[2][j] + (gamma - 1.0)*(qnew[2][j] - 0.5*qnew[1][j]*qnew[1][j]/qnew[0][j]))/qnew[0][j];
  137. 		}
  138.  
  139. 	}
  140.  
  141. 	//ファイル出力//
  142. 	for(j=0;j<100;j++){
  143. 		fprintf(fp,"%.4f ", count[j]);
  144. 	}
  145. 	fprintf(fp,"\n");
  146. 	for(j=0;j<100;j++){
  147. 		fprintf(fp,"%.4f ", rho[j]);
  148. 	}
  149. 	fprintf(fp,"\n");
  150. 	for(j=0;j<100;j++){
  151. 		fprintf(fp,"%.4f ", u[j]);
  152. 	}
  153. 	fprintf(fp,"\n");
  154. 	for(j=0;j<100;j++){
  155. 		fprintf(fp,"%.4f ", p[j]);
  156. 	}
  157. 	fclose(fp);
  158.  
  159. 	return 0;
  160. }
Time limit exceeded #stdin #stdout 5s 2808KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/5ccG50
language:
C++ (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!