miniWeather_openacc.cpp 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646
  1. //////////////////////////////////////////////////////////////////////////////////////////
  2. // miniWeather
  3. // Author: Matt Norman <normanmr@ornl.gov> , Oak Ridge National Laboratory
  4. // This code simulates dry, stratified, compressible, non-hydrostatic fluid flows
  5. // For documentation, please see the attached documentation in the "documentation" folder
  6. //////////////////////////////////////////////////////////////////////////////////////////
  7. /*
  8. ** Copyright (c) 2018, National Center for Computational Sciences, Oak Ridge National Laboratory. All rights reserved.
  9. **
  10. ** Portions Copyright (c) 2020, NVIDIA Corporation. All rights reserved.
  11. */
  12. #include <stdlib.h>
  13. #include <math.h>
  14. #include <stdio.h>
  15. #include <nvtx3/nvToolsExt.h>
  16. const double pi = 3.14159265358979323846264338327; //Pi
  17. const double grav = 9.8; //Gravitational acceleration (m / s^2)
  18. const double cp = 1004.; //Specific heat of dry air at constant pressure
  19. const double rd = 287.; //Dry air constant for equation of state (P=rho*rd*T)
  20. const double p0 = 1.e5; //Standard pressure at the surface in Pascals
  21. const double C0 = 27.5629410929725921310572974482; //Constant to translate potential temperature into pressure (P=C0*(rho*theta)**gamma)
  22. const double gamm = 1.40027894002789400278940027894; //gamma=cp/Rd , have to call this gamm because "gamma" is taken (I hate C so much)
  23. //Define domain and stability-related constants
  24. const double xlen = 2.e4; //Length of the domain in the x-direction (meters)
  25. const double zlen = 1.e4; //Length of the domain in the z-direction (meters)
  26. const double hv_beta = 0.25; //How strong to diffuse the solution: hv_beta \in [0:1]
  27. const double cfl = 1.50; //"Courant, Friedrichs, Lewy" number (for numerical stability)
  28. const double max_speed = 450; //Assumed maximum wave speed during the simulation (speed of sound + speed of wind) (meter / sec)
  29. const int hs = 2; //"Halo" size: number of cells needed for a full "stencil" of information for reconstruction
  30. const int sten_size = 4; //Size of the stencil used for interpolation
  31. //Parameters for indexing and flags
  32. const int NUM_VARS = 4; //Number of fluid state variables
  33. const int ID_DENS = 0; //index for density ("rho")
  34. const int ID_UMOM = 1; //index for momentum in the x-direction ("rho * u")
  35. const int ID_WMOM = 2; //index for momentum in the z-direction ("rho * w")
  36. const int ID_RHOT = 3; //index for density * potential temperature ("rho * theta")
  37. const int DIR_X = 1; //Integer constant to express that this operation is in the x-direction
  38. const int DIR_Z = 2; //Integer constant to express that this operation is in the z-direction
  39. const int nqpoints = 3;
  40. double qpoints[] = {0.112701665379258311482073460022E0, 0.500000000000000000000000000000E0, 0.887298334620741688517926539980E0};
  41. double qweights[] = {0.277777777777777777777777777779E0, 0.444444444444444444444444444444E0, 0.277777777777777777777777777779E0};
  42. ///////////////////////////////////////////////////////////////////////////////////////
  43. // Variables that are initialized but remain static over the course of the simulation
  44. ///////////////////////////////////////////////////////////////////////////////////////
  45. double sim_time; //total simulation time in seconds
  46. double output_freq; //frequency to perform output in seconds
  47. double dt; //Model time step (seconds)
  48. int nx, nz; //Number of local grid cells in the x- and z- dimensions
  49. double dx, dz; //Grid space length in x- and z-dimension (meters)
  50. int nx_glob, nz_glob; //Number of total grid cells in the x- and z- dimensions
  51. int i_beg, k_beg; //beginning index in the x- and z-directions
  52. int nranks, myrank; //my rank id
  53. int left_rank, right_rank; //Rank IDs that exist to my left and right in the global domain
  54. double *hy_dens_cell; //hydrostatic density (vert cell avgs). Dimensions: (1-hs:nz+hs)
  55. double *hy_dens_theta_cell; //hydrostatic rho*t (vert cell avgs). Dimensions: (1-hs:nz+hs)
  56. double *hy_dens_int; //hydrostatic density (vert cell interf). Dimensions: (1:nz+1)
  57. double *hy_dens_theta_int; //hydrostatic rho*t (vert cell interf). Dimensions: (1:nz+1)
  58. double *hy_pressure_int; //hydrostatic press (vert cell interf). Dimensions: (1:nz+1)
  59. ///////////////////////////////////////////////////////////////////////////////////////
  60. // Variables that are dynamics over the course of the simulation
  61. ///////////////////////////////////////////////////////////////////////////////////////
  62. double etime; //Elapsed model time
  63. double output_counter; //Helps determine when it's time to do output
  64. //Runtime variable arrays
  65. double *state; //Fluid state. Dimensions: (1-hs:nx+hs,1-hs:nz+hs,NUM_VARS)
  66. double *state_tmp; //Fluid state. Dimensions: (1-hs:nx+hs,1-hs:nz+hs,NUM_VARS)
  67. double *flux; //Cell interface fluxes. Dimensions: (nx+1,nz+1,NUM_VARS)
  68. double *tend; //Fluid state tendencies. Dimensions: (nx,nz,NUM_VARS)
  69. int num_out = 0; //The number of outputs performed so far
  70. int direction_switch = 1;
  71. //How is this not in the standard?!
  72. double dmin(double a, double b)
  73. {
  74. if (a < b)
  75. {
  76. return a;
  77. }
  78. else
  79. {
  80. return b;
  81. }
  82. };
  83. //Declaring the functions defined after "main"
  84. void init();
  85. void finalize();
  86. void injection(double x, double z, double &r, double &u, double &w, double &t, double &hr, double &ht);
  87. void hydro_const_theta(double z, double &r, double &t);
  88. void output(double *state, double etime);
  89. void ncwrap(int ierr, int line);
  90. void perform_timestep(double *state, double *state_tmp, double *flux, double *tend, double dt);
  91. void semi_discrete_step(double *state_init, double *state_forcing, double *state_out, double dt, int dir, double *flux, double *tend);
  92. void compute_tendencies_x(double *state, double *flux, double *tend);
  93. void compute_tendencies_z(double *state, double *flux, double *tend);
  94. void set_halo_values_x(double *state);
  95. void set_halo_values_z(double *state);
  96. ///////////////////////////////////////////////////////////////////////////////////////
  97. // THE MAIN PROGRAM STARTS HERE
  98. ///////////////////////////////////////////////////////////////////////////////////////
  99. int main(int argc, char **argv)
  100. {
  101. ///////////////////////////////////////////////////////////////////////////////////////
  102. // BEGIN USER-CONFIGURABLE PARAMETERS
  103. ///////////////////////////////////////////////////////////////////////////////////////
  104. //The x-direction length is twice as long as the z-direction length
  105. //So, you'll want to have nx_glob be twice as large as nz_glob
  106. nx_glob = 40; //Number of total cells in the x-dirction
  107. nz_glob = 20; //Number of total cells in the z-dirction
  108. sim_time = 1000; //How many seconds to run the simulation
  109. output_freq = 100; //How frequently to output data to file (in seconds)
  110. ///////////////////////////////////////////////////////////////////////////////////////
  111. // END USER-CONFIGURABLE PARAMETERS
  112. ///////////////////////////////////////////////////////////////////////////////////////
  113. if (argc == 4)
  114. {
  115. printf("The arguments supplied are %s %s %s\n", argv[1], argv[2], argv[3]);
  116. nx_glob = atoi(argv[1]);
  117. nz_glob = atoi(argv[2]);
  118. sim_time = atoi(argv[3]);
  119. }
  120. else
  121. {
  122. printf("Using default values ...\n");
  123. }
  124. nvtxRangePushA("Total");
  125. init();
  126. //Output the initial state
  127. //output(state, etime);
  128. ////////////////////////////////////////////////////
  129. // MAIN TIME STEP LOOP
  130. ////////////////////////////////////////////////////
  131. nvtxRangePushA("while");
  132. while (etime < sim_time)
  133. {
  134. //If the time step leads to exceeding the simulation time, shorten it for the last step
  135. if (etime + dt > sim_time)
  136. {
  137. dt = sim_time - etime;
  138. }
  139. //Perform a single time step
  140. nvtxRangePushA("perform_timestep");
  141. perform_timestep(state, state_tmp, flux, tend, dt);
  142. nvtxRangePop();
  143. //Inform the user
  144. printf("Elapsed Time: %lf / %lf\n", etime, sim_time);
  145. //Update the elapsed time and output counter
  146. etime = etime + dt;
  147. output_counter = output_counter + dt;
  148. //If it's time for output, reset the counter, and do output
  149. if (output_counter >= output_freq)
  150. {
  151. output_counter = output_counter - output_freq;
  152. //output(state, etime);
  153. }
  154. }
  155. nvtxRangePop();
  156. finalize();
  157. nvtxRangePop();
  158. }
  159. //Performs a single dimensionally split time step using a simple low-storate three-stage Runge-Kutta time integrator
  160. //The dimensional splitting is a second-order-accurate alternating Strang splitting in which the
  161. //order of directions is alternated each time step.
  162. //The Runge-Kutta method used here is defined as follows:
  163. // q* = q[n] + dt/3 * rhs(q[n])
  164. // q** = q[n] + dt/2 * rhs(q* )
  165. // q[n+1] = q[n] + dt/1 * rhs(q** )
  166. void perform_timestep(double *state, double *state_tmp, double *flux, double *tend, double dt)
  167. {
  168. if (direction_switch)
  169. {
  170. //x-direction first
  171. semi_discrete_step(state, state, state_tmp, dt / 3, DIR_X, flux, tend);
  172. semi_discrete_step(state, state_tmp, state_tmp, dt / 2, DIR_X, flux, tend);
  173. semi_discrete_step(state, state_tmp, state, dt / 1, DIR_X, flux, tend);
  174. //z-direction second
  175. semi_discrete_step(state, state, state_tmp, dt / 3, DIR_Z, flux, tend);
  176. semi_discrete_step(state, state_tmp, state_tmp, dt / 2, DIR_Z, flux, tend);
  177. semi_discrete_step(state, state_tmp, state, dt / 1, DIR_Z, flux, tend);
  178. }
  179. else
  180. {
  181. //z-direction second
  182. semi_discrete_step(state, state, state_tmp, dt / 3, DIR_Z, flux, tend);
  183. semi_discrete_step(state, state_tmp, state_tmp, dt / 2, DIR_Z, flux, tend);
  184. semi_discrete_step(state, state_tmp, state, dt / 1, DIR_Z, flux, tend);
  185. //x-direction first
  186. semi_discrete_step(state, state, state_tmp, dt / 3, DIR_X, flux, tend);
  187. semi_discrete_step(state, state_tmp, state_tmp, dt / 2, DIR_X, flux, tend);
  188. semi_discrete_step(state, state_tmp, state, dt / 1, DIR_X, flux, tend);
  189. }
  190. if (direction_switch)
  191. {
  192. direction_switch = 0;
  193. }
  194. else
  195. {
  196. direction_switch = 1;
  197. }
  198. }
  199. //Perform a single semi-discretized step in time with the form:
  200. //state_out = state_init + dt * rhs(state_forcing)
  201. //Meaning the step starts from state_init, computes the rhs using state_forcing, and stores the result in state_out
  202. void semi_discrete_step(double *state_init, double *state_forcing, double *state_out, double dt, int dir, double *flux, double *tend)
  203. {
  204. int i, k, ll, inds, indt;
  205. if (dir == DIR_X)
  206. {
  207. //Set the halo values in the x-direction
  208. set_halo_values_x(state_forcing);
  209. //Compute the time tendencies for the fluid state in the x-direction
  210. compute_tendencies_x(state_forcing, flux, tend);
  211. }
  212. else if (dir == DIR_Z)
  213. {
  214. //Set the halo values in the z-direction
  215. set_halo_values_z(state_forcing);
  216. //Compute the time tendencies for the fluid state in the z-direction
  217. compute_tendencies_z(state_forcing, flux, tend);
  218. }
  219. /////////////////////////////////////////////////
  220. // TODO: THREAD ME
  221. /////////////////////////////////////////////////
  222. //Apply the tendencies to the fluid state
  223. #pragma acc parallel loop private(inds, indt) default(present)
  224. for (ll = 0; ll < NUM_VARS; ll++)
  225. {
  226. for (k = 0; k < nz; k++)
  227. {
  228. for (i = 0; i < nx; i++)
  229. {
  230. inds = ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + i + hs;
  231. indt = ll * nz * nx + k * nx + i;
  232. state_out[inds] = state_init[inds] + dt * tend[indt];
  233. }
  234. }
  235. }
  236. }
  237. //Compute the time tendencies of the fluid state using forcing in the x-direction
  238. //First, compute the flux vector at each cell interface in the x-direction (including hyperviscosity)
  239. //Then, compute the tendencies using those fluxes
  240. void compute_tendencies_x(double *state, double *flux, double *tend)
  241. {
  242. int i, k, ll, s, inds, indf1, indf2, indt;
  243. double r, u, w, t, p, stencil[4], d3_vals[NUM_VARS], vals[NUM_VARS], hv_coef;
  244. //Compute the hyperviscosity coeficient
  245. hv_coef = -hv_beta * dx / (16 * dt);
  246. /////////////////////////////////////////////////
  247. // TODO: THREAD ME
  248. /////////////////////////////////////////////////
  249. //Compute fluxes in the x-direction for each cell
  250. #pragma acc parallel loop private(ll, s, inds, stencil, vals, d3_vals, r, u, w, t, p)
  251. for (k = 0; k < nz; k++)
  252. {
  253. for (i = 0; i < nx + 1; i++)
  254. {
  255. //Use fourth-order interpolation from four cell averages to compute the value at the interface in question
  256. for (ll = 0; ll < NUM_VARS; ll++)
  257. {
  258. for (s = 0; s < sten_size; s++)
  259. {
  260. inds = ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + i + s;
  261. stencil[s] = state[inds];
  262. }
  263. //Fourth-order-accurate interpolation of the state
  264. vals[ll] = -stencil[0] / 12 + 7 * stencil[1] / 12 + 7 * stencil[2] / 12 - stencil[3] / 12;
  265. //First-order-accurate interpolation of the third spatial derivative of the state (for artificial viscosity)
  266. d3_vals[ll] = -stencil[0] + 3 * stencil[1] - 3 * stencil[2] + stencil[3];
  267. }
  268. //Compute density, u-wind, w-wind, potential temperature, and pressure (r,u,w,t,p respectively)
  269. r = vals[ID_DENS] + hy_dens_cell[k + hs];
  270. u = vals[ID_UMOM] / r;
  271. w = vals[ID_WMOM] / r;
  272. t = (vals[ID_RHOT] + hy_dens_theta_cell[k + hs]) / r;
  273. p = C0 * pow((r * t), gamm);
  274. //Compute the flux vector
  275. flux[ID_DENS * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * u - hv_coef * d3_vals[ID_DENS];
  276. flux[ID_UMOM * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * u * u + p - hv_coef * d3_vals[ID_UMOM];
  277. flux[ID_WMOM * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * u * w - hv_coef * d3_vals[ID_WMOM];
  278. flux[ID_RHOT * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * u * t - hv_coef * d3_vals[ID_RHOT];
  279. }
  280. }
  281. /////////////////////////////////////////////////
  282. // TODO: THREAD ME
  283. /////////////////////////////////////////////////
  284. //Use the fluxes to compute tendencies for each cell
  285. #pragma acc parallel loop private(indt, indf1, indf2)
  286. for (ll = 0; ll < NUM_VARS; ll++)
  287. {
  288. for (k = 0; k < nz; k++)
  289. {
  290. for (i = 0; i < nx; i++)
  291. {
  292. indt = ll * nz * nx + k * nx + i;
  293. indf1 = ll * (nz + 1) * (nx + 1) + k * (nx + 1) + i;
  294. indf2 = ll * (nz + 1) * (nx + 1) + k * (nx + 1) + i + 1;
  295. tend[indt] = -(flux[indf2] - flux[indf1]) / dx;
  296. }
  297. }
  298. }
  299. }
  300. //Compute the time tendencies of the fluid state using forcing in the z-direction
  301. //First, compute the flux vector at each cell interface in the z-direction (including hyperviscosity)
  302. //Then, compute the tendencies using those fluxes
  303. void compute_tendencies_z(double *state, double *flux, double *tend)
  304. {
  305. int i, k, ll, s, inds, indf1, indf2, indt;
  306. double r, u, w, t, p, stencil[4], d3_vals[NUM_VARS], vals[NUM_VARS], hv_coef;
  307. //Compute the hyperviscosity coeficient
  308. hv_coef = -hv_beta * dx / (16 * dt);
  309. /////////////////////////////////////////////////
  310. // TODO: THREAD ME
  311. /////////////////////////////////////////////////
  312. //Compute fluxes in the x-direction for each cell
  313. #pragma acc parallel loop private(ll, s, inds, stencil, vals, d3_vals, r, u, w, t, p)
  314. for (k = 0; k < nz + 1; k++)
  315. {
  316. for (i = 0; i < nx; i++)
  317. {
  318. //Use fourth-order interpolation from four cell averages to compute the value at the interface in question
  319. for (ll = 0; ll < NUM_VARS; ll++)
  320. {
  321. for (s = 0; s < sten_size; s++)
  322. {
  323. inds = ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + s) * (nx + 2 * hs) + i + hs;
  324. stencil[s] = state[inds];
  325. }
  326. //Fourth-order-accurate interpolation of the state
  327. vals[ll] = -stencil[0] / 12 + 7 * stencil[1] / 12 + 7 * stencil[2] / 12 - stencil[3] / 12;
  328. //First-order-accurate interpolation of the third spatial derivative of the state
  329. d3_vals[ll] = -stencil[0] + 3 * stencil[1] - 3 * stencil[2] + stencil[3];
  330. }
  331. //Compute density, u-wind, w-wind, potential temperature, and pressure (r,u,w,t,p respectively)
  332. r = vals[ID_DENS] + hy_dens_int[k];
  333. u = vals[ID_UMOM] / r;
  334. w = vals[ID_WMOM] / r;
  335. t = (vals[ID_RHOT] + hy_dens_theta_int[k]) / r;
  336. p = C0 * pow((r * t), gamm) - hy_pressure_int[k];
  337. //Compute the flux vector with hyperviscosity
  338. flux[ID_DENS * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * w - hv_coef * d3_vals[ID_DENS];
  339. flux[ID_UMOM * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * w * u - hv_coef * d3_vals[ID_UMOM];
  340. flux[ID_WMOM * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * w * w + p - hv_coef * d3_vals[ID_WMOM];
  341. flux[ID_RHOT * (nz + 1) * (nx + 1) + k * (nx + 1) + i] = r * w * t - hv_coef * d3_vals[ID_RHOT];
  342. }
  343. }
  344. /////////////////////////////////////////////////
  345. // TODO: THREAD ME
  346. /////////////////////////////////////////////////
  347. //Use the fluxes to compute tendencies for each cell
  348. #pragma acc parallel loop private(indt, indf1, indf2)
  349. for (ll = 0; ll < NUM_VARS; ll++)
  350. {
  351. for (k = 0; k < nz; k++)
  352. {
  353. for (i = 0; i < nx; i++)
  354. {
  355. indt = ll * nz * nx + k * nx + i;
  356. indf1 = ll * (nz + 1) * (nx + 1) + (k) * (nx + 1) + i;
  357. indf2 = ll * (nz + 1) * (nx + 1) + (k + 1) * (nx + 1) + i;
  358. tend[indt] = -(flux[indf2] - flux[indf1]) / dz;
  359. if (ll == ID_WMOM)
  360. {
  361. inds = ID_DENS * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + i + hs;
  362. tend[indt] = tend[indt] - state[inds] * grav;
  363. }
  364. }
  365. }
  366. }
  367. }
  368. void set_halo_values_x(double *state)
  369. {
  370. int k, ll, ind_r, ind_u, ind_t, i;
  371. double z;
  372. #pragma acc parallel loop
  373. for (ll = 0; ll < NUM_VARS; ll++)
  374. {
  375. for (k = 0; k < nz; k++)
  376. {
  377. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + 0] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + nx + hs - 2];
  378. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + 1] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + nx + hs - 1];
  379. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + nx + hs] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + hs];
  380. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + nx + hs + 1] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + hs + 1];
  381. }
  382. }
  383. ////////////////////////////////////////////////////
  384. if (myrank == 0)
  385. {
  386. for (k = 0; k < nz; k++)
  387. {
  388. for (i = 0; i < hs; i++)
  389. {
  390. z = (k_beg + k + 0.5) * dz;
  391. if (abs(z - 3 * zlen / 4) <= zlen / 16)
  392. {
  393. ind_r = ID_DENS * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + i;
  394. ind_u = ID_UMOM * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + i;
  395. ind_t = ID_RHOT * (nz + 2 * hs) * (nx + 2 * hs) + (k + hs) * (nx + 2 * hs) + i;
  396. state[ind_u] = (state[ind_r] + hy_dens_cell[k + hs]) * 50.;
  397. state[ind_t] = (state[ind_r] + hy_dens_cell[k + hs]) * 298. - hy_dens_theta_cell[k + hs];
  398. }
  399. }
  400. }
  401. }
  402. }
  403. //Set this task's halo values in the z-direction.
  404. //decomposition in the vertical direction.
  405. void set_halo_values_z(double *state)
  406. {
  407. int i, ll;
  408. const double mnt_width = xlen / 8;
  409. double x, xloc, mnt_deriv;
  410. /////////////////////////////////////////////////
  411. // TODO: THREAD ME
  412. /////////////////////////////////////////////////
  413. #pragma acc parallel loop private(x, xloc, mnt_deriv)
  414. for (ll = 0; ll < NUM_VARS; ll++)
  415. {
  416. for (i = 0; i < nx + 2 * hs; i++)
  417. {
  418. if (ll == ID_WMOM)
  419. {
  420. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (0) * (nx + 2 * hs) + i] = 0.;
  421. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (1) * (nx + 2 * hs) + i] = 0.;
  422. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (nz + hs) * (nx + 2 * hs) + i] = 0.;
  423. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (nz + hs + 1) * (nx + 2 * hs) + i] = 0.;
  424. }
  425. else
  426. {
  427. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (0) * (nx + 2 * hs) + i] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (hs) * (nx + 2 * hs) + i];
  428. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (1) * (nx + 2 * hs) + i] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (hs) * (nx + 2 * hs) + i];
  429. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (nz + hs) * (nx + 2 * hs) + i] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (nz + hs - 1) * (nx + 2 * hs) + i];
  430. state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (nz + hs + 1) * (nx + 2 * hs) + i] = state[ll * (nz + 2 * hs) * (nx + 2 * hs) + (nz + hs - 1) * (nx + 2 * hs) + i];
  431. }
  432. }
  433. }
  434. }
  435. void init()
  436. {
  437. int i, k, ii, kk, ll, inds, i_end;
  438. double x, z, r, u, w, t, hr, ht, nper;
  439. //Set the cell grid size
  440. dx = xlen / nx_glob;
  441. dz = zlen / nz_glob;
  442. nranks = 1;
  443. myrank = 0;
  444. // For simpler version, replace i_beg = 0, nx = nx_glob, left_rank = 0, right_rank = 0;
  445. nper = ((double)nx_glob) / nranks;
  446. i_beg = round(nper * (myrank));
  447. i_end = round(nper * ((myrank) + 1)) - 1;
  448. nx = i_end - i_beg + 1;
  449. left_rank = myrank - 1;
  450. if (left_rank == -1)
  451. left_rank = nranks - 1;
  452. right_rank = myrank + 1;
  453. if (right_rank == nranks)
  454. right_rank = 0;
  455. ////////////////////////////////////////////////////////////////////////////////
  456. ////////////////////////////////////////////////////////////////////////////////
  457. // YOU DON'T NEED TO ALTER ANYTHING BELOW THIS POINT IN THE CODE
  458. ////////////////////////////////////////////////////////////////////////////////
  459. ////////////////////////////////////////////////////////////////////////////////
  460. k_beg = 0;
  461. nz = nz_glob;
  462. //Allocate the model data
  463. state = (double *)malloc((nx + 2 * hs) * (nz + 2 * hs) * NUM_VARS * sizeof(double));
  464. state_tmp = (double *)malloc((nx + 2 * hs) * (nz + 2 * hs) * NUM_VARS * sizeof(double));
  465. flux = (double *)malloc((nx + 1) * (nz + 1) * NUM_VARS * sizeof(double));
  466. tend = (double *)malloc(nx * nz * NUM_VARS * sizeof(double));
  467. hy_dens_cell = (double *)malloc((nz + 2 * hs) * sizeof(double));
  468. hy_dens_theta_cell = (double *)malloc((nz + 2 * hs) * sizeof(double));
  469. hy_dens_int = (double *)malloc((nz + 1) * sizeof(double));
  470. hy_dens_theta_int = (double *)malloc((nz + 1) * sizeof(double));
  471. hy_pressure_int = (double *)malloc((nz + 1) * sizeof(double));
  472. //Define the maximum stable time step based on an assumed maximum wind speed
  473. dt = dmin(dx, dz) / max_speed * cfl;
  474. //Set initial elapsed model time and output_counter to zero
  475. etime = 0.;
  476. output_counter = 0.;
  477. // Display grid information
  478. printf("nx_glob, nz_glob: %d %d\n", nx_glob, nz_glob);
  479. printf("dx,dz: %lf %lf\n", dx, dz);
  480. printf("dt: %lf\n", dt);
  481. //////////////////////////////////////////////////////////////////////////
  482. // Initialize the cell-averaged fluid state via Gauss-Legendre quadrature
  483. //////////////////////////////////////////////////////////////////////////
  484. for (k = 0; k < nz + 2 * hs; k++)
  485. {
  486. for (i = 0; i < nx + 2 * hs; i++)
  487. {
  488. //Initialize the state to zero
  489. for (ll = 0; ll < NUM_VARS; ll++)
  490. {
  491. inds = ll * (nz + 2 * hs) * (nx + 2 * hs) + k * (nx + 2 * hs) + i;
  492. state[inds] = 0.;
  493. }
  494. //Use Gauss-Legendre quadrature to initialize a hydrostatic balance + temperature perturbation
  495. for (kk = 0; kk < nqpoints; kk++)
  496. {
  497. for (ii = 0; ii < nqpoints; ii++)
  498. {
  499. //Compute the x,z location within the global domain based on cell and quadrature index
  500. x = (i_beg + i - hs + 0.5) * dx + (qpoints[ii] - 0.5) * dx;
  501. z = (k_beg + k - hs + 0.5) * dz + (qpoints[kk] - 0.5) * dz;
  502. //Set the fluid state based on the user's specification (default is injection in this example)
  503. injection(x, z, r, u, w, t, hr, ht);
  504. //Store into the fluid state array
  505. inds = ID_DENS * (nz + 2 * hs) * (nx + 2 * hs) + k * (nx + 2 * hs) + i;
  506. state[inds] = state[inds] + r * qweights[ii] * qweights[kk];
  507. inds = ID_UMOM * (nz + 2 * hs) * (nx + 2 * hs) + k * (nx + 2 * hs) + i;
  508. state[inds] = state[inds] + (r + hr) * u * qweights[ii] * qweights[kk];
  509. inds = ID_WMOM * (nz + 2 * hs) * (nx + 2 * hs) + k * (nx + 2 * hs) + i;
  510. state[inds] = state[inds] + (r + hr) * w * qweights[ii] * qweights[kk];
  511. inds = ID_RHOT * (nz + 2 * hs) * (nx + 2 * hs) + k * (nx + 2 * hs) + i;
  512. state[inds] = state[inds] + ((r + hr) * (t + ht) - hr * ht) * qweights[ii] * qweights[kk];
  513. }
  514. }
  515. for (ll = 0; ll < NUM_VARS; ll++)
  516. {
  517. inds = ll * (nz + 2 * hs) * (nx + 2 * hs) + k * (nx + 2 * hs) + i;
  518. state_tmp[inds] = state[inds];
  519. }
  520. }
  521. }
  522. //Compute the hydrostatic background state over vertical cell averages
  523. for (k = 0; k < nz + 2 * hs; k++)
  524. {
  525. hy_dens_cell[k] = 0.;
  526. hy_dens_theta_cell[k] = 0.;
  527. for (kk = 0; kk < nqpoints; kk++)
  528. {
  529. z = (k_beg + k - hs + 0.5) * dz;
  530. //Set the fluid state based on the user's specification (default is injection in this example)
  531. injection(0., z, r, u, w, t, hr, ht);
  532. hy_dens_cell[k] = hy_dens_cell[k] + hr * qweights[kk];
  533. hy_dens_theta_cell[k] = hy_dens_theta_cell[k] + hr * ht * qweights[kk];
  534. }
  535. }
  536. //Compute the hydrostatic background state at vertical cell interfaces
  537. for (k = 0; k < nz + 1; k++)
  538. {
  539. z = (k_beg + k) * dz;
  540. //Set the fluid state based on the user's specification (default is injection in this example)
  541. injection(0., z, r, u, w, t, hr, ht);
  542. hy_dens_int[k] = hr;
  543. hy_dens_theta_int[k] = hr * ht;
  544. hy_pressure_int[k] = C0 * pow((hr * ht), gamm);
  545. }
  546. }
  547. //This test case is initially balanced but injects fast, cold air from the left boundary near the model top
  548. //x and z are input coordinates at which to sample
  549. //r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
  550. //hr and ht are output background hydrostatic density and potential temperature at that location
  551. void injection(double x, double z, double &r, double &u, double &w, double &t, double &hr, double &ht)
  552. {
  553. hydro_const_theta(z, hr, ht);
  554. r = 0.;
  555. t = 0.;
  556. u = 0.;
  557. w = 0.;
  558. }
  559. //Establish hydrstatic balance using constant potential temperature (thermally neutral atmosphere)
  560. //z is the input coordinate
  561. //r and t are the output background hydrostatic density and potential temperature
  562. void hydro_const_theta(double z, double &r, double &t)
  563. {
  564. const double theta0 = 300.; //Background potential temperature
  565. const double exner0 = 1.; //Surface-level Exner pressure
  566. double p, exner, rt;
  567. //Establish hydrostatic balance first using Exner pressure
  568. t = theta0; //Potential Temperature at z
  569. exner = exner0 - grav * z / (cp * theta0); //Exner pressure at z
  570. p = p0 * pow(exner, (cp / rd)); //Pressure at z
  571. rt = pow((p / C0), (1. / gamm)); //rho*theta at z
  572. r = rt / t; //Density at z
  573. }
  574. void finalize()
  575. {
  576. free(state);
  577. free(state_tmp);
  578. free(flux);
  579. free(tend);
  580. free(hy_dens_cell);
  581. free(hy_dens_theta_cell);
  582. free(hy_dens_int);
  583. free(hy_dens_theta_int);
  584. free(hy_pressure_int);
  585. }