API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Feb 21st, 2018
61
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 32.58 KB | None | 0 0
raw download clone embed print report
  1. #include <chrono>
  2. #include <iostream>
  3. #include <math.h>
  4. #include <numeric>
  5.  
  6. // Global constants in the equations are
  7. double rkold[3];
  8. int nx1;
  9. double rinv8;
  10. double rinv9;
  11. double Minf;
  12. double rinv1;
  13. double rinv4;
  14. double rinv5;
  15. double Pr;
  16. double rinv12;
  17. double rinv13;
  18. double rinv10;
  19. double rinv11;
  20. double rknew[3];
  21. double rc6;
  22. double rc7;
  23. double rc0;
  24. double rc2;
  25. double rc3;
  26. int nx0;
  27. double deltai1;
  28. double deltai0;
  29. double Re;
  30. double deltat;
  31. double gama;
  32. int itercount;
  33.  
  34. // main program start
  35. int main(int argc, char **argv) {
  36.  
  37. // Initialising global constants
  38. nx0 = 257;
  39. if (argc > 1)
  40. nx0 = atoi(argv[1]);
  41. nx1 = 257;
  42. if (argc > 2)
  43. nx1 = atoi(argv[2]);
  44. gama = 1.40000000000000;
  45. Pr = 0.710000000000000;
  46. Re = 1600;
  47. deltat = 0.000846250000000000;
  48. Minf = 0.100000000000000;
  49. rc6 = 5.0 / 2.0;
  50. rc7 = 4.0 / 3.0;
  51. rc0 = 1.0 / 2.0;
  52. rc2 = 1.0 / 12.0;
  53. rc3 = 2.0 / 3.0;
  54. rkold[0] = 1.0 / 4.0;
  55. rkold[1] = 3.0 / 20.0;
  56. rkold[2] = 3.0 / 5.0;
  57. rknew[0] = 2.0 / 3.0;
  58. rknew[1] = 5.0 / 12.0;
  59. rknew[2] = 3.0 / 5.0;
  60. rinv12 = pow(Minf, -2);
  61. rinv13 = 1.0 / (gama * pow(Minf, 2));
  62. rinv10 = 1.0 / Pr;
  63. rinv11 = 1.0 / (gama - 1);
  64. rinv9 = 1.0 / Re;
  65. deltai1 = (1.0 / (nx1 - 1.0)) * M_PI;
  66. deltai0 = (1.0 / (nx0 - 1.0)) * M_PI;
  67. rinv8 = pow(deltai1, -2);
  68. rinv1 = 1.0 / deltai1;
  69. rinv4 = 1.0 / deltai0;
  70. rinv5 = pow(deltai0, -2);
  71. itercount = 20;
  72.  
  73. std::cout << "Running on a " << nx0 + 4 << "x" << nx1 + 4 << " mesh for "
  74. << itercount << " iterations\n";
  75.  
  76. // Allocating mesh
  77. double *rho = new double[(nx0 + 4) * (nx1 + 4)];
  78. double *rhou0 = new double[(nx0 + 4) * (nx1 + 4)];
  79. double *rhou1 = new double[(nx0 + 4) * (nx1 + 4)];
  80. double *rhoE = new double[(nx0 + 4) * (nx1 + 4)];
  81. double *rho_old = new double[(nx0 + 4) * (nx1 + 4)];
  82. double *rhou0_old = new double[(nx0 + 4) * (nx1 + 4)];
  83. double *rhou1_old = new double[(nx0 + 4) * (nx1 + 4)];
  84. double *rhoE_old = new double[(nx0 + 4) * (nx1 + 4)];
  85. double *T = new double[(nx0 + 4) * (nx1 + 4)];
  86. double *u0 = new double[(nx0 + 4) * (nx1 + 4)];
  87. double *u1 = new double[(nx0 + 4) * (nx1 + 4)];
  88. double *p = new double[(nx0 + 4) * (nx1 + 4)];
  89. double *wk0 = new double[(nx0 + 4) * (nx1 + 4)];
  90. double *wk1 = new double[(nx0 + 4) * (nx1 + 4)];
  91. double *wk2 = new double[(nx0 + 4) * (nx1 + 4)];
  92. double *wk3 = new double[(nx0 + 4) * (nx1 + 4)];
  93.  
  94. // Initialisation
  95. //writing dataset rho with (j, i)
  96. //writing dataset rhou0 with (i, j)
  97. //writing dataset rhou1 with (i, j)
  98. //writing dataset rhoE with (i, j)
  99. for (int j = 0; j < nx1 + 4; j++) {
  100. for (int i = 0; i < nx0 + 4; i++) {
  101. double x = deltai0 * (i - 2);
  102. double y = deltai1 * (j - 2);
  103. double u = sin(x) * cos(y);
  104. double v = -cos(x) * sin(y);
  105. double p = 1.0 * rinv13 + 0.25 * (sin(2.0 * x) + sin(2.0 * y));
  106. double r = gama * pow(Minf, 2) * p;
  107. rho[(j + 0) * (nx0 + 4) + (i + 0)] = r;
  108. rhou0[(j + 0) * (nx0 + 4) + (i + 0)] = r * u;
  109. rhou1[(j + 0) * (nx0 + 4) + (i + 0)] = r * v;
  110. rhoE[(j + 0) * (nx0 + 4) + (i + 0)] =
  111. rinv11 * p + 0.5 * r * (pow(u, 2) + pow(v, 2));
  112. }
  113. }
  114.  
  115. // Apply boundary conditions
  116. // Left
  117. //writing dataset rho with (i-1, j), (i-2, j)
  118. //reading dataset rho with (i+1, j), (i+2, j)
  119. //writing dataset rhou0 with (i-1, j), (i-2, j)
  120. //reading dataset rhou0 with (i+1, j), (i+2, j)
  121. //writing dataset rhou1 with (i-1, j), (i-2, j)
  122. //reading dataset rhou1 with (i+1, j), (i+2, j)
  123. //writing dataset rhouE with (i-1, j), (i-2, j)
  124. //reading dataset rhouE with (i+1, j), (i+2, j)
  125.  
  126. for (int j = 0; j < nx1 + 4; j++) {
  127. for (int i = 2; i < 3; i++) {
  128. rho[(j + 0) * (nx0 + 4) + (i - 1)] = rho[(j + 0) * (nx0 + 4) + (i + 1)];
  129. rho[(j + 0) * (nx0 + 4) + (i - 2)] = rho[(j + 0) * (nx0 + 4) + (i + 2)];
  130. rhou0[(j + 0) * (nx0 + 4) + (i - 1)] =
  131. rhou0[(j + 0) * (nx0 + 4) + (i + 1)];
  132. rhou0[(j + 0) * (nx0 + 4) + (i - 2)] =
  133. rhou0[(j + 0) * (nx0 + 4) + (i + 2)];
  134. rhou1[(j + 0) * (nx0 + 4) + (i - 1)] =
  135. rhou1[(j + 0) * (nx0 + 4) + (i + 1)];
  136. rhou1[(j + 0) * (nx0 + 4) + (i - 2)] =
  137. rhou1[(j + 0) * (nx0 + 4) + (i + 2)];
  138. rhoE[(j + 0) * (nx0 + 4) + (i - 1)] = rhoE[(j + 0) * (nx0 + 4) + (i + 1)];
  139. rhoE[(j + 0) * (nx0 + 4) + (i - 2)] = rhoE[(j + 0) * (nx0 + 4) + (i + 2)];
  140. }
  141. }
  142.  
  143. // Right
  144. //writing dataset rho with (i+1, j), (i+2, j)
  145. //reading dataset rho with (i-1, j), (i-2, j)
  146. //writing dataset rhou0 with (i+1, j), (i+2, j)
  147. //reading dataset rhou0 with (i-1, j), (i-2, j)
  148. for (int j = 0; j < nx1 + 4; j++) {
  149. for (int i = nx0 + 1; i < nx0 + 2; i++) {
  150. rho[(j + 0) * (nx0 + 4) + (i + 1)] = rho[(j + 0) * (nx0 + 4) + (i - 1)];
  151. rho[(j + 0) * (nx0 + 4) + (i + 2)] = rho[(j + 0) * (nx0 + 4) + (i - 2)];
  152. rhou0[(j + 0) * (nx0 + 4) + (i + 1)] =
  153. rhou0[(j + 0) * (nx0 + 4) + (i - 1)];
  154. rhou0[(j + 0) * (nx0 + 4) + (i + 2)] =
  155. rhou0[(j + 0) * (nx0 + 4) + (i - 2)];
  156. rhou1[(j + 0) * (nx0 + 4) + (i + 1)] =
  157. rhou1[(j + 0) * (nx0 + 4) + (i - 1)];
  158. rhou1[(j + 0) * (nx0 + 4) + (i + 2)] =
  159. rhou1[(j + 0) * (nx0 + 4) + (i - 2)];
  160. rhoE[(j + 0) * (nx0 + 4) + (i + 1)] = rhoE[(j + 0) * (nx0 + 4) + (i - 1)];
  161. rhoE[(j + 0) * (nx0 + 4) + (i + 2)] = rhoE[(j + 0) * (nx0 + 4) + (i - 2)];
  162. }
  163. }
  164.  
  165. // Top
  166. for (int j = 2; j < 3; j++) {
  167. for (int i = 0; i < nx0 + 4; i++) {
  168. rho[(j - 1) * (nx0 + 4) + (i + 0)] = rho[(j + 1) * (nx0 + 4) + (i + 0)];
  169. rho[(j - 2) * (nx0 + 4) + (i + 0)] = rho[(j + 2) * (nx0 + 4) + (i + 0)];
  170. rhou0[(j - 1) * (nx0 + 4) + (i + 0)] =
  171. rhou0[(j + 1) * (nx0 + 4) + (i + 0)];
  172. rhou0[(j - 2) * (nx0 + 4) + (i + 0)] =
  173. rhou0[(j + 2) * (nx0 + 4) + (i + 0)];
  174. rhou1[(j - 1) * (nx0 + 4) + (i + 0)] =
  175. rhou1[(j + 1) * (nx0 + 4) + (i + 0)];
  176. rhou1[(j - 2) * (nx0 + 4) + (i + 0)] =
  177. rhou1[(j + 2) * (nx0 + 4) + (i + 0)];
  178. rhoE[(j - 1) * (nx0 + 4) + (i + 0)] = rhoE[(j + 1) * (nx0 + 4) + (i + 0)];
  179. rhoE[(j - 2) * (nx0 + 4) + (i + 0)] = rhoE[(j + 2) * (nx0 + 4) + (i + 0)];
  180. }
  181. }
  182.  
  183. // Bottom
  184. for (int j = nx1 + 1; j < nx1 + 2; j++) {
  185. for (int i = 0; i < nx0 + 4; i++) {
  186. rho[(j + 1) * (nx0 + 4) + (i + 0)] = rho[(j - 1) * (nx0 + 4) + (i + 0)];
  187. rho[(j + 2) * (nx0 + 4) + (i + 0)] = rho[(j - 2) * (nx0 + 4) + (i + 0)];
  188. rhou0[(j + 1) * (nx0 + 4) + (i + 0)] =
  189. rhou0[(j - 1) * (nx0 + 4) + (i + 0)];
  190. rhou0[(j + 2) * (nx0 + 4) + (i + 0)] =
  191. rhou0[(j - 2) * (nx0 + 4) + (i + 0)];
  192. rhou1[(j + 1) * (nx0 + 4) + (i + 0)] =
  193. rhou1[(j - 1) * (nx0 + 4) + (i + 0)];
  194. rhou1[(j + 2) * (nx0 + 4) + (i + 0)] =
  195. rhou1[(j - 2) * (nx0 + 4) + (i + 0)];
  196. rhoE[(j + 1) * (nx0 + 4) + (i + 0)] = rhoE[(j - 1) * (nx0 + 4) + (i + 0)];
  197. rhoE[(j + 2) * (nx0 + 4) + (i + 0)] = rhoE[(j - 2) * (nx0 + 4) + (i + 0)];
  198. }
  199. }
  200.  
  201. // Record start time
  202. auto start = std::chrono::high_resolution_clock::now();
  203.  
  204. // Main time iteration loop
  205. for (int iteration = 0; iteration < itercount; iteration++) {
  206.  
  207. // Save equations
  208. for (int j = 0; j < nx1 + 4; j++) {
  209. for (int i = 0; i < nx0 + 4; i++) {
  210. rho_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  211. rho[(j + 0) * (nx0 + 4) + (i + 0)];
  212. rhou0_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  213. rhou0[(j + 0) * (nx0 + 4) + (i + 0)];
  214. rhou1_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  215. rhou1[(j + 0) * (nx0 + 4) + (i + 0)];
  216. rhoE_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  217. rhoE[(j + 0) * (nx0 + 4) + (i + 0)];
  218. }
  219. }
  220.  
  221. // Runge-Kutta time-stepper
  222. for (int stage = 0; stage < 3; stage++) {
  223.  
  224. // Grouped Formula Evaluation
  225. for (int j = 0; j < nx1 + 4; j++) {
  226. for (int i = 0; i < nx0 + 4; i++) {
  227. T[(j + 0) * (nx0 + 4) + (i + 0)] =
  228. gama * (gama - 1) *
  229. ((-rc0 * pow(rhou0[(j + 0) * (nx0 + 4) + (i + 0)], 2) -
  230. rc0 * pow(rhou1[(j + 0) * (nx0 + 4) + (i + 0)], 2)) /
  231. rho[(j + 0) * (nx0 + 4) + (i + 0)] +
  232. rhoE[(j + 0) * (nx0 + 4) + (i + 0)]) *
  233. pow(Minf, 2) / rho[(j + 0) * (nx0 + 4) + (i + 0)];
  234. p[(j + 0) * (nx0 + 4) + (i + 0)] =
  235. (gama - 1) *
  236. ((-rc0 * pow(rhou0[(j + 0) * (nx0 + 4) + (i + 0)], 2) -
  237. rc0 * pow(rhou1[(j + 0) * (nx0 + 4) + (i + 0)], 2)) /
  238. rho[(j + 0) * (nx0 + 4) + (i + 0)] +
  239. rhoE[(j + 0) * (nx0 + 4) + (i + 0)]);
  240. u1[(j + 0) * (nx0 + 4) + (i + 0)] =
  241. rhou1[(j + 0) * (nx0 + 4) + (i + 0)] /
  242. rho[(j + 0) * (nx0 + 4) + (i + 0)];
  243. u0[(j + 0) * (nx0 + 4) + (i + 0)] =
  244. rhou0[(j + 0) * (nx0 + 4) + (i + 0)] /
  245. rho[(j + 0) * (nx0 + 4) + (i + 0)];
  246. }
  247. }
  248.  
  249. // Residual of equation
  250. for (int j = 2; j < nx1 + 2; j++) {
  251. for (int i = 2; i < nx0 + 2; i++) {
  252. double temp_eval0 =
  253. rinv1 * ((rc2)*rhoE[(j - 2) * (nx0 + 4) + (i + 0)] *
  254. u1[(j - 2) * (nx0 + 4) + (i + 0)] -
  255. rc3 * rhoE[(j - 1) * (nx0 + 4) + (i + 0)] *
  256. u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  257. (rc3)*rhoE[(j + 1) * (nx0 + 4) + (i + 0)] *
  258. u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  259. rc2 * rhoE[(j + 2) * (nx0 + 4) + (i + 0)] *
  260. u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  261. double temp_eval1 =
  262. rinv4 * ((rc2)*rhoE[(j + 0) * (nx0 + 4) + (i - 2)] -
  263. rc3 * rhoE[(j + 0) * (nx0 + 4) + (i - 1)] +
  264. (rc3)*rhoE[(j + 0) * (nx0 + 4) + (i + 1)] -
  265. rc2 * rhoE[(j + 0) * (nx0 + 4) + (i + 2)]);
  266. double temp_eval2 =
  267. rinv4 * ((rc2)*rhoE[(j + 0) * (nx0 + 4) + (i - 2)] *
  268. u0[(j + 0) * (nx0 + 4) + (i - 2)] -
  269. rc3 * rhoE[(j + 0) * (nx0 + 4) + (i - 1)] *
  270. u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  271. (rc3)*rhoE[(j + 0) * (nx0 + 4) + (i + 1)] *
  272. u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  273. rc2 * rhoE[(j + 0) * (nx0 + 4) + (i + 2)] *
  274. u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  275. double temp_eval3 = rinv1 * ((rc2)*u1[(j - 2) * (nx0 + 4) + (i + 0)] -
  276. rc3 * u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  277. (rc3)*u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  278. rc2 * u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  279. double temp_eval4 = rinv5 * (-rc6 * T[(j + 0) * (nx0 + 4) + (i + 0)] -
  280. rc2 * T[(j + 0) * (nx0 + 4) + (i - 2)] +
  281. (rc7)*T[(j + 0) * (nx0 + 4) + (i - 1)] +
  282. (rc7)*T[(j + 0) * (nx0 + 4) + (i + 1)] -
  283. rc2 * T[(j + 0) * (nx0 + 4) + (i + 2)]);
  284. double temp_eval5 =
  285. rinv1 * ((rc2)*rhou1[(j - 2) * (nx0 + 4) + (i + 0)] -
  286. rc3 * rhou1[(j - 1) * (nx0 + 4) + (i + 0)] +
  287. (rc3)*rhou1[(j + 1) * (nx0 + 4) + (i + 0)] -
  288. rc2 * rhou1[(j + 2) * (nx0 + 4) + (i + 0)]);
  289. double temp_eval6 =
  290. rinv5 * (-rc6 * u0[(j + 0) * (nx0 + 4) + (i + 0)] -
  291. rc2 * u0[(j + 0) * (nx0 + 4) + (i - 2)] +
  292. (rc7)*u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  293. (rc7)*u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  294. rc2 * u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  295. double temp_eval7 =
  296. rinv1 * ((rc2)*rhou1[(j - 2) * (nx0 + 4) + (i + 0)] *
  297. u1[(j - 2) * (nx0 + 4) + (i + 0)] -
  298. rc3 * rhou1[(j - 1) * (nx0 + 4) + (i + 0)] *
  299. u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  300. (rc3)*rhou1[(j + 1) * (nx0 + 4) + (i + 0)] *
  301. u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  302. rc2 * rhou1[(j + 2) * (nx0 + 4) + (i + 0)] *
  303. u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  304. double temp_eval8 =
  305. rinv1 * ((rc2)*rho[(j - 2) * (nx0 + 4) + (i + 0)] -
  306. rc3 * rho[(j - 1) * (nx0 + 4) + (i + 0)] +
  307. (rc3)*rho[(j + 1) * (nx0 + 4) + (i + 0)] -
  308. rc2 * rho[(j + 2) * (nx0 + 4) + (i + 0)]);
  309. double temp_eval9 =
  310. rinv4 * ((rc2)*rhou0[(j + 0) * (nx0 + 4) + (i - 2)] -
  311. rc3 * rhou0[(j + 0) * (nx0 + 4) + (i - 1)] +
  312. (rc3)*rhou0[(j + 0) * (nx0 + 4) + (i + 1)] -
  313. rc2 * rhou0[(j + 0) * (nx0 + 4) + (i + 2)]);
  314. double temp_eval10 =
  315. rinv1 * ((rc2)*u0[(j - 2) * (nx0 + 4) + (i + 0)] -
  316. rc3 * u0[(j - 1) * (nx0 + 4) + (i + 0)] +
  317. (rc3)*u0[(j + 1) * (nx0 + 4) + (i + 0)] -
  318. rc2 * u0[(j + 2) * (nx0 + 4) + (i + 0)]);
  319. double temp_eval11 = rinv1 * ((rc2)*p[(j - 2) * (nx0 + 4) + (i + 0)] -
  320. rc3 * p[(j - 1) * (nx0 + 4) + (i + 0)] +
  321. (rc3)*p[(j + 1) * (nx0 + 4) + (i + 0)] -
  322. rc2 * p[(j + 2) * (nx0 + 4) + (i + 0)]);
  323. double temp_eval12 =
  324. rinv4 * ((rc2)*u1[(j + 0) * (nx0 + 4) + (i - 2)] -
  325. rc3 * u1[(j + 0) * (nx0 + 4) + (i - 1)] +
  326. (rc3)*u1[(j + 0) * (nx0 + 4) + (i + 1)] -
  327. rc2 * u1[(j + 0) * (nx0 + 4) + (i + 2)]);
  328. double temp_eval13 =
  329. rinv4 * ((rc2)*rhou1[(j + 0) * (nx0 + 4) + (i - 2)] *
  330. u0[(j + 0) * (nx0 + 4) + (i - 2)] -
  331. rc3 * rhou1[(j + 0) * (nx0 + 4) + (i - 1)] *
  332. u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  333. (rc3)*rhou1[(j + 0) * (nx0 + 4) + (i + 1)] *
  334. u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  335. rc2 * rhou1[(j + 0) * (nx0 + 4) + (i + 2)] *
  336. u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  337. double temp_eval14 =
  338. rinv4 * ((rc2)*rho[(j + 0) * (nx0 + 4) + (i - 2)] *
  339. u0[(j + 0) * (nx0 + 4) + (i - 2)] -
  340. rc3 * rho[(j + 0) * (nx0 + 4) + (i - 1)] *
  341. u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  342. (rc3)*rho[(j + 0) * (nx0 + 4) + (i + 1)] *
  343. u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  344. rc2 * rho[(j + 0) * (nx0 + 4) + (i + 2)] *
  345. u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  346. double temp_eval15 =
  347. rinv4 * ((rc2)*rho[(j + 0) * (nx0 + 4) + (i - 2)] -
  348. rc3 * rho[(j + 0) * (nx0 + 4) + (i - 1)] +
  349. (rc3)*rho[(j + 0) * (nx0 + 4) + (i + 1)] -
  350. rc2 * rho[(j + 0) * (nx0 + 4) + (i + 2)]);
  351. double temp_eval16 =
  352. rinv1 * ((rc2)*rhou0[(j - 2) * (nx0 + 4) + (i + 0)] -
  353. rc3 * rhou0[(j - 1) * (nx0 + 4) + (i + 0)] +
  354. (rc3)*rhou0[(j + 1) * (nx0 + 4) + (i + 0)] -
  355. rc2 * rhou0[(j + 2) * (nx0 + 4) + (i + 0)]);
  356. double temp_eval17 =
  357. rinv5 * (-rc6 * u1[(j + 0) * (nx0 + 4) + (i + 0)] -
  358. rc2 * u1[(j + 0) * (nx0 + 4) + (i - 2)] +
  359. (rc7)*u1[(j + 0) * (nx0 + 4) + (i - 1)] +
  360. (rc7)*u1[(j + 0) * (nx0 + 4) + (i + 1)] -
  361. rc2 * u1[(j + 0) * (nx0 + 4) + (i + 2)]);
  362. double temp_eval18 =
  363. rinv4 * ((rc2)*rhou0[(j + 0) * (nx0 + 4) + (i - 2)] *
  364. u0[(j + 0) * (nx0 + 4) + (i - 2)] -
  365. rc3 * rhou0[(j + 0) * (nx0 + 4) + (i - 1)] *
  366. u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  367. (rc3)*rhou0[(j + 0) * (nx0 + 4) + (i + 1)] *
  368. u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  369. rc2 * rhou0[(j + 0) * (nx0 + 4) + (i + 2)] *
  370. u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  371. double temp_eval19 =
  372. rinv1 * ((rc2)*rhou0[(j - 2) * (nx0 + 4) + (i + 0)] *
  373. u1[(j - 2) * (nx0 + 4) + (i + 0)] -
  374. rc3 * rhou0[(j - 1) * (nx0 + 4) + (i + 0)] *
  375. u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  376. (rc3)*rhou0[(j + 1) * (nx0 + 4) + (i + 0)] *
  377. u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  378. rc2 * rhou0[(j + 2) * (nx0 + 4) + (i + 0)] *
  379. u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  380. double temp_eval20 =
  381. rinv8 * (-rc6 * u1[(j + 0) * (nx0 + 4) + (i + 0)] -
  382. rc2 * u1[(j - 2) * (nx0 + 4) + (i + 0)] +
  383. (rc7)*u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  384. (rc7)*u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  385. rc2 * u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  386. double temp_eval21 = rinv4 * ((rc2)*p[(j + 0) * (nx0 + 4) + (i - 2)] -
  387. rc3 * p[(j + 0) * (nx0 + 4) + (i - 1)] +
  388. (rc3)*p[(j + 0) * (nx0 + 4) + (i + 1)] -
  389. rc2 * p[(j + 0) * (nx0 + 4) + (i + 2)]);
  390. double temp_eval22 =
  391. rinv8 * (-rc6 * u0[(j + 0) * (nx0 + 4) + (i + 0)] -
  392. rc2 * u0[(j - 2) * (nx0 + 4) + (i + 0)] +
  393. (rc7)*u0[(j - 1) * (nx0 + 4) + (i + 0)] +
  394. (rc7)*u0[(j + 1) * (nx0 + 4) + (i + 0)] -
  395. rc2 * u0[(j + 2) * (nx0 + 4) + (i + 0)]);
  396. double temp_eval23 =
  397. rinv1 * ((rc2)*rhoE[(j - 2) * (nx0 + 4) + (i + 0)] -
  398. rc3 * rhoE[(j - 1) * (nx0 + 4) + (i + 0)] +
  399. (rc3)*rhoE[(j + 1) * (nx0 + 4) + (i + 0)] -
  400. rc2 * rhoE[(j + 2) * (nx0 + 4) + (i + 0)]);
  401. double temp_eval24 =
  402. rinv8 * (-rc6 * T[(j + 0) * (nx0 + 4) + (i + 0)] -
  403. rc2 * T[(j - 2) * (nx0 + 4) + (i + 0)] +
  404. (rc7)*T[(j - 1) * (nx0 + 4) + (i + 0)] +
  405. (rc7)*T[(j + 1) * (nx0 + 4) + (i + 0)] -
  406. rc2 * T[(j + 2) * (nx0 + 4) + (i + 0)]);
  407. double temp_eval25 =
  408. rinv4 * ((rc2)*rhou1[(j + 0) * (nx0 + 4) + (i - 2)] -
  409. rc3 * rhou1[(j + 0) * (nx0 + 4) + (i - 1)] +
  410. (rc3)*rhou1[(j + 0) * (nx0 + 4) + (i + 1)] -
  411. rc2 * rhou1[(j + 0) * (nx0 + 4) + (i + 2)]);
  412. double temp_eval26 = rinv4 * ((rc2)*p[(j + 0) * (nx0 + 4) + (i - 2)] *
  413. u0[(j + 0) * (nx0 + 4) + (i - 2)] -
  414. rc3 * p[(j + 0) * (nx0 + 4) + (i - 1)] *
  415. u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  416. (rc3)*p[(j + 0) * (nx0 + 4) + (i + 1)] *
  417. u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  418. rc2 * p[(j + 0) * (nx0 + 4) + (i + 2)] *
  419. u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  420. double temp_eval27 =
  421. rinv4 * ((rc2)*u0[(j + 0) * (nx0 + 4) + (i - 2)] -
  422. rc3 * u0[(j + 0) * (nx0 + 4) + (i - 1)] +
  423. (rc3)*u0[(j + 0) * (nx0 + 4) + (i + 1)] -
  424. rc2 * u0[(j + 0) * (nx0 + 4) + (i + 2)]);
  425. double temp_eval28 = rinv1 * ((rc2)*p[(j - 2) * (nx0 + 4) + (i + 0)] *
  426. u1[(j - 2) * (nx0 + 4) + (i + 0)] -
  427. rc3 * p[(j - 1) * (nx0 + 4) + (i + 0)] *
  428. u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  429. (rc3)*p[(j + 1) * (nx0 + 4) + (i + 0)] *
  430. u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  431. rc2 * p[(j + 2) * (nx0 + 4) + (i + 0)] *
  432. u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  433. double temp_eval29 =
  434. rinv1 * ((rc2)*rho[(j - 2) * (nx0 + 4) + (i + 0)] *
  435. u1[(j - 2) * (nx0 + 4) + (i + 0)] -
  436. rc3 * rho[(j - 1) * (nx0 + 4) + (i + 0)] *
  437. u1[(j - 1) * (nx0 + 4) + (i + 0)] +
  438. (rc3)*rho[(j + 1) * (nx0 + 4) + (i + 0)] *
  439. u1[(j + 1) * (nx0 + 4) + (i + 0)] -
  440. rc2 * rho[(j + 2) * (nx0 + 4) + (i + 0)] *
  441. u1[(j + 2) * (nx0 + 4) + (i + 0)]);
  442. double temp_eval30 =
  443. rinv1 * ((rc2)*rinv4 * ((rc2)*u0[(j - 2) * (nx0 + 4) + (i - 2)] -
  444. rc3 * u0[(j - 2) * (nx0 + 4) + (i - 1)] +
  445. (rc3)*u0[(j - 2) * (nx0 + 4) + (i + 1)] -
  446. rc2 * u0[(j - 2) * (nx0 + 4) + (i + 2)]) -
  447. rc3 * rinv4 * ((rc2)*u0[(j - 1) * (nx0 + 4) + (i - 2)] -
  448. rc3 * u0[(j - 1) * (nx0 + 4) + (i - 1)] +
  449. (rc3)*u0[(j - 1) * (nx0 + 4) + (i + 1)] -
  450. rc2 * u0[(j - 1) * (nx0 + 4) + (i + 2)]) +
  451. (rc3)*rinv4 * ((rc2)*u0[(j + 1) * (nx0 + 4) + (i - 2)] -
  452. rc3 * u0[(j + 1) * (nx0 + 4) + (i - 1)] +
  453. (rc3)*u0[(j + 1) * (nx0 + 4) + (i + 1)] -
  454. rc2 * u0[(j + 1) * (nx0 + 4) + (i + 2)]) -
  455. rc2 * rinv4 * ((rc2)*u0[(j + 2) * (nx0 + 4) + (i - 2)] -
  456. rc3 * u0[(j + 2) * (nx0 + 4) + (i - 1)] +
  457. (rc3)*u0[(j + 2) * (nx0 + 4) + (i + 1)] -
  458. rc2 * u0[(j + 2) * (nx0 + 4) + (i + 2)]));
  459. double temp_eval31 =
  460. rinv1 * ((rc2)*rinv4 * ((rc2)*u1[(j - 2) * (nx0 + 4) + (i - 2)] -
  461. rc3 * u1[(j - 2) * (nx0 + 4) + (i - 1)] +
  462. (rc3)*u1[(j - 2) * (nx0 + 4) + (i + 1)] -
  463. rc2 * u1[(j - 2) * (nx0 + 4) + (i + 2)]) -
  464. rc3 * rinv4 * ((rc2)*u1[(j - 1) * (nx0 + 4) + (i - 2)] -
  465. rc3 * u1[(j - 1) * (nx0 + 4) + (i - 1)] +
  466. (rc3)*u1[(j - 1) * (nx0 + 4) + (i + 1)] -
  467. rc2 * u1[(j - 1) * (nx0 + 4) + (i + 2)]) +
  468. (rc3)*rinv4 * ((rc2)*u1[(j + 1) * (nx0 + 4) + (i - 2)] -
  469. rc3 * u1[(j + 1) * (nx0 + 4) + (i - 1)] +
  470. (rc3)*u1[(j + 1) * (nx0 + 4) + (i + 1)] -
  471. rc2 * u1[(j + 1) * (nx0 + 4) + (i + 2)]) -
  472. rc2 * rinv4 * ((rc2)*u1[(j + 2) * (nx0 + 4) + (i - 2)] -
  473. rc3 * u1[(j + 2) * (nx0 + 4) + (i - 1)] +
  474. (rc3)*u1[(j + 2) * (nx0 + 4) + (i + 1)] -
  475. rc2 * u1[(j + 2) * (nx0 + 4) + (i + 2)]));
  476. wk0[(j + 0) * (nx0 + 4) + (i + 0)] =
  477. -0.5 * temp_eval14 -
  478. 0.5 * temp_eval15 * u0[(j + 0) * (nx0 + 4) + (i + 0)] -
  479. 0.5 * temp_eval29 -
  480. 0.5 * temp_eval8 * u1[(j + 0) * (nx0 + 4) + (i + 0)] -
  481. 0.5 * (temp_eval27 + temp_eval3) *
  482. rho[(j + 0) * (nx0 + 4) + (i + 0)];
  483. wk1[(j + 0) * (nx0 + 4) + (i + 0)] =
  484. -0.5 * temp_eval16 * u1[(j + 0) * (nx0 + 4) + (i + 0)] -
  485. 0.5 * temp_eval18 - 0.5 * temp_eval19 - temp_eval21 -
  486. 0.5 * temp_eval9 * u0[(j + 0) * (nx0 + 4) + (i + 0)] +
  487. rinv9 * (temp_eval22 + temp_eval31) +
  488. rinv9 * (-rc3 * temp_eval31 + (rc7)*temp_eval6) -
  489. 0.5 * (temp_eval27 + temp_eval3) *
  490. rhou0[(j + 0) * (nx0 + 4) + (i + 0)];
  491. wk2[(j + 0) * (nx0 + 4) + (i + 0)] =
  492. -temp_eval11 - 0.5 * temp_eval13 -
  493. 0.5 * temp_eval25 * u0[(j + 0) * (nx0 + 4) + (i + 0)] -
  494. 0.5 * temp_eval5 * u1[(j + 0) * (nx0 + 4) + (i + 0)] -
  495. 0.5 * temp_eval7 + rinv9 * (temp_eval17 + temp_eval30) +
  496. rinv9 * ((rc7)*temp_eval20 - rc3 * temp_eval30) -
  497. 0.5 * (temp_eval27 + temp_eval3) *
  498. rhou1[(j + 0) * (nx0 + 4) + (i + 0)];
  499. wk3[(j + 0) * (nx0 + 4) + (i + 0)] =
  500. -0.5 * temp_eval0 -
  501. 0.5 * temp_eval1 * u0[(j + 0) * (nx0 + 4) + (i + 0)] +
  502. temp_eval10 * rinv9 * (temp_eval10 + temp_eval12) +
  503. temp_eval12 * rinv9 * (temp_eval10 + temp_eval12) -
  504. 0.5 * temp_eval2 -
  505. 0.5 * temp_eval23 * u1[(j + 0) * (nx0 + 4) + (i + 0)] +
  506. temp_eval24 * rinv10 * rinv11 * rinv12 * rinv9 - temp_eval26 +
  507. temp_eval27 * rinv9 * ((rc7)*temp_eval27 - rc3 * temp_eval3) -
  508. temp_eval28 +
  509. temp_eval3 * rinv9 * (-rc3 * temp_eval27 + (rc7)*temp_eval3) +
  510. temp_eval4 * rinv10 * rinv11 * rinv12 * rinv9 +
  511. rinv9 * (temp_eval17 + temp_eval30) *
  512. u1[(j + 0) * (nx0 + 4) + (i + 0)] +
  513. rinv9 * ((rc7)*temp_eval20 - rc3 * temp_eval30) *
  514. u1[(j + 0) * (nx0 + 4) + (i + 0)] +
  515. rinv9 * (temp_eval22 + temp_eval31) *
  516. u0[(j + 0) * (nx0 + 4) + (i + 0)] +
  517. rinv9 * (-rc3 * temp_eval31 + (rc7)*temp_eval6) *
  518. u0[(j + 0) * (nx0 + 4) + (i + 0)] -
  519. 0.5 * (temp_eval27 + temp_eval3) *
  520. rhoE[(j + 0) * (nx0 + 4) + (i + 0)];
  521. }
  522. }
  523.  
  524. // RK new (subloop) update
  525. for (int j = 0; j < nx1 + 4; j++) {
  526. for (int i = 0; i < nx0 + 4; i++) {
  527. rho[(j + 0) * (nx0 + 4) + (i + 0)] =
  528. deltat * rknew[0] * wk0[(j + 0) * (nx0 + 4) + (i + 0)] +
  529. rho_old[(j + 0) * (nx0 + 4) + (i + 0)];
  530. rhou0[(j + 0) * (nx0 + 4) + (i + 0)] =
  531. deltat * rknew[0] * wk1[(j + 0) * (nx0 + 4) + (i + 0)] +
  532. rhou0_old[(j + 0) * (nx0 + 4) + (i + 0)];
  533. rhou1[(j + 0) * (nx0 + 4) + (i + 0)] =
  534. deltat * rknew[0] * wk2[(j + 0) * (nx0 + 4) + (i + 0)] +
  535. rhou1_old[(j + 0) * (nx0 + 4) + (i + 0)];
  536. rhoE[(j + 0) * (nx0 + 4) + (i + 0)] =
  537. deltat * rknew[0] * wk3[(j + 0) * (nx0 + 4) + (i + 0)] +
  538. rhoE_old[(j + 0) * (nx0 + 4) + (i + 0)];
  539. }
  540. }
  541.  
  542. // RK old update
  543. for (int j = 0; j < nx1 + 4; j++) {
  544. for (int i = 0; i < nx0 + 4; i++) {
  545. rho_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  546. deltat * rkold[0] * wk0[(j + 0) * (nx0 + 4) + (i + 0)] +
  547. rho_old[(j + 0) * (nx0 + 4) + (i + 0)];
  548. rhou0_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  549. deltat * rkold[0] * wk1[(j + 0) * (nx0 + 4) + (i + 0)] +
  550. rhou0_old[(j + 0) * (nx0 + 4) + (i + 0)];
  551. rhou1_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  552. deltat * rkold[0] * wk2[(j + 0) * (nx0 + 4) + (i + 0)] +
  553. rhou1_old[(j + 0) * (nx0 + 4) + (i + 0)];
  554. rhoE_old[(j + 0) * (nx0 + 4) + (i + 0)] =
  555. deltat * rkold[0] * wk3[(j + 0) * (nx0 + 4) + (i + 0)] +
  556. rhoE_old[(j + 0) * (nx0 + 4) + (i + 0)];
  557. }
  558. }
  559.  
  560. // Apply boundary conditions
  561.  
  562. // Left
  563. for (int j = 0; j < nx1 + 4; j++) {
  564. for (int i = 2; i < 3; i++) {
  565. rho[(j + 0) * (nx0 + 4) + (i - 1)] =
  566. rho[(j + 0) * (nx0 + 4) + (i + 1)];
  567. rho[(j + 0) * (nx0 + 4) + (i - 2)] =
  568. rho[(j + 0) * (nx0 + 4) + (i + 2)];
  569. rhou0[(j + 0) * (nx0 + 4) + (i - 1)] =
  570. rhou0[(j + 0) * (nx0 + 4) + (i + 1)];
  571. rhou0[(j + 0) * (nx0 + 4) + (i - 2)] =
  572. rhou0[(j + 0) * (nx0 + 4) + (i + 2)];
  573. rhou1[(j + 0) * (nx0 + 4) + (i - 1)] =
  574. rhou1[(j + 0) * (nx0 + 4) + (i + 1)];
  575. rhou1[(j + 0) * (nx0 + 4) + (i - 2)] =
  576. rhou1[(j + 0) * (nx0 + 4) + (i + 2)];
  577. rhoE[(j + 0) * (nx0 + 4) + (i - 1)] =
  578. rhoE[(j + 0) * (nx0 + 4) + (i + 1)];
  579. rhoE[(j + 0) * (nx0 + 4) + (i - 2)] =
  580. rhoE[(j + 0) * (nx0 + 4) + (i + 2)];
  581. }
  582. }
  583.  
  584. // Right
  585. for (int j = 0; j < nx1 + 4; j++) {
  586. for (int i = nx0 + 1; i < nx0 + 2; i++) {
  587. rho[(j + 0) * (nx0 + 4) + (i + 1)] =
  588. rho[(j + 0) * (nx0 + 4) + (i - 1)];
  589. rho[(j + 0) * (nx0 + 4) + (i + 2)] =
  590. rho[(j + 0) * (nx0 + 4) + (i - 2)];
  591. rhou0[(j + 0) * (nx0 + 4) + (i + 1)] =
  592. rhou0[(j + 0) * (nx0 + 4) + (i - 1)];
  593. rhou0[(j + 0) * (nx0 + 4) + (i + 2)] =
  594. rhou0[(j + 0) * (nx0 + 4) + (i - 2)];
  595. rhou1[(j + 0) * (nx0 + 4) + (i + 1)] =
  596. rhou1[(j + 0) * (nx0 + 4) + (i - 1)];
  597. rhou1[(j + 0) * (nx0 + 4) + (i + 2)] =
  598. rhou1[(j + 0) * (nx0 + 4) + (i - 2)];
  599. rhoE[(j + 0) * (nx0 + 4) + (i + 1)] =
  600. rhoE[(j + 0) * (nx0 + 4) + (i - 1)];
  601. rhoE[(j + 0) * (nx0 + 4) + (i + 2)] =
  602. rhoE[(j + 0) * (nx0 + 4) + (i - 2)];
  603. }
  604. }
  605.  
  606. // Top
  607. for (int j = 2; j < 3; j++) {
  608. for (int i = 0; i < nx0 + 4; i++) {
  609. rho[(j - 1) * (nx0 + 4) + (i + 0)] =
  610. rho[(j + 1) * (nx0 + 4) + (i + 0)];
  611. rho[(j - 2) * (nx0 + 4) + (i + 0)] =
  612. rho[(j + 2) * (nx0 + 4) + (i + 0)];
  613. rhou0[(j - 1) * (nx0 + 4) + (i + 0)] =
  614. rhou0[(j + 1) * (nx0 + 4) + (i + 0)];
  615. rhou0[(j - 2) * (nx0 + 4) + (i + 0)] =
  616. rhou0[(j + 2) * (nx0 + 4) + (i + 0)];
  617. rhou1[(j - 1) * (nx0 + 4) + (i + 0)] =
  618. rhou1[(j + 1) * (nx0 + 4) + (i + 0)];
  619. rhou1[(j - 2) * (nx0 + 4) + (i + 0)] =
  620. rhou1[(j + 2) * (nx0 + 4) + (i + 0)];
  621. rhoE[(j - 1) * (nx0 + 4) + (i + 0)] =
  622. rhoE[(j + 1) * (nx0 + 4) + (i + 0)];
  623. rhoE[(j - 2) * (nx0 + 4) + (i + 0)] =
  624. rhoE[(j + 2) * (nx0 + 4) + (i + 0)];
  625. }
  626. }
  627.  
  628. // Bottom
  629. for (int j = nx1 + 1; j < nx1 + 2; j++) {
  630. for (int i = 0; i < nx0 + 4; i++) {
  631. rho[(j + 1) * (nx0 + 4) + (i + 0)] =
  632. rho[(j - 1) * (nx0 + 4) + (i + 0)];
  633. rho[(j + 2) * (nx0 + 4) + (i + 0)] =
  634. rho[(j - 2) * (nx0 + 4) + (i + 0)];
  635. rhou0[(j + 1) * (nx0 + 4) + (i + 0)] =
  636. rhou0[(j - 1) * (nx0 + 4) + (i + 0)];
  637. rhou0[(j + 2) * (nx0 + 4) + (i + 0)] =
  638. rhou0[(j - 2) * (nx0 + 4) + (i + 0)];
  639. rhou1[(j + 1) * (nx0 + 4) + (i + 0)] =
  640. rhou1[(j - 1) * (nx0 + 4) + (i + 0)];
  641. rhou1[(j + 2) * (nx0 + 4) + (i + 0)] =
  642. rhou1[(j - 2) * (nx0 + 4) + (i + 0)];
  643. rhoE[(j + 1) * (nx0 + 4) + (i + 0)] =
  644. rhoE[(j - 1) * (nx0 + 4) + (i + 0)];
  645. rhoE[(j + 2) * (nx0 + 4) + (i + 0)] =
  646. rhoE[(j - 2) * (nx0 + 4) + (i + 0)];
  647. }
  648. }
  649. } // End of stage loop
  650.  
  651. double sum = 0.0;
  652. double sum2 = 0.0;
  653. for (int j = 0; j < nx1 + 4; j++) {
  654. for (int i = 0; i < nx0 + 4; i++) {
  655. sum += rho[j * (nx0 + 4) + i] * rho[j * (nx0 + 4) + i];
  656. sum2 += p[j * (nx0 + 4) + i] * p[j * (nx0 + 4) + i];
  657. }
  658. }
  659. std::cout << "Checksums: " << sqrt(sum) << " " << sqrt(sum2) << "\n";
  660.  
  661. } // End of time loop
  662.  
  663. // Record end time
  664. auto end = std::chrono::high_resolution_clock::now();
  665. std::chrono::duration<double> diff = end - start;
  666.  
  667. std::cout << "\nTimings are:\n";
  668. std::cout << "-----------------------------------------\n";
  669. // TODO: per-loop statistics come here
  670. std::cout << "Total Wall time " << diff.count() << " seconds\n";
  671.  
  672. delete[] rho;
  673. delete[] rhou0;
  674. delete[] rhou1;
  675. delete[] rhoE;
  676. delete[] rho_old;
  677. delete[] rhou0_old;
  678. delete[] rhou1_old;
  679. delete[] rhoE_old;
  680. delete[] T;
  681. delete[] u0;
  682. delete[] u1;
  683. delete[] p;
  684. delete[] wk0;
  685. delete[] wk1;
  686. delete[] wk2;
  687. delete[] wk3;
  688. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!