Canoe
Comprehensive Atmosphere N' Ocean Engine
netcdf.cpp
Go to the documentation of this file.
1 // C/C++
2 #include <cstdio>
3 #include <cstdlib>
4 #include <iomanip>
5 #include <iostream>
6 #include <sstream>
7 #include <stdexcept>
8 #include <string>
9 
10 // canoe
11 #include <configure.hpp>
12 #include <impl.hpp>
13 
14 // athena
15 #include <athena/athena.hpp>
16 #include <athena/coordinates/coordinates.hpp>
17 #include <athena/mesh/mesh.hpp>
18 #include <athena/outputs/user_outputs.hpp>
19 
20 // harp
21 #include <harp/radiation.hpp>
22 #include <harp/radiation_band.hpp>
23 
24 // output
25 #include "output_utils.hpp"
26 
27 // Only proceed if NETCDF output enabled
28 #ifdef NETCDFOUTPUT
29 
30 // External library headers
31 #include <netcdf.h>
32 
33 //----------------------------------------------------------------------------------------
34 // NetcdfOutput constructor
35 // destructor - not needed for this derived class
36 
37 NetcdfOutput::NetcdfOutput(OutputParameters oparams) : OutputType(oparams) {}
38 
39 //----------------------------------------------------------------------------------------
41 // \brief Cycles over all MeshBlocks and writes OutputData in NETCDF format,
42 // one
43 // MeshBlock per file
44 
45 void NetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
46  auto pmeta = MetadataTable::GetInstance();
47 
48  // Loop over MeshBlocks
49  for (int b = 0; b < pm->nblocal; ++b) {
50  MeshBlock *pmb = pm->my_blocks(b);
51  auto prad = pmb->pimpl->prad;
52 
53  // set start/end array indices depending on whether ghost zones are included
54  out_is = pmb->is;
55  out_ie = pmb->ie;
56  out_js = pmb->js;
57  out_je = pmb->je;
58  out_ks = pmb->ks;
59  out_ke = pmb->ke;
60 
61  // FIXME: include_ghost zones probably doesn't work with grids other than
62  // CCC
63  if (output_params.include_ghost_zones) {
64  out_is -= NGHOST;
65  out_ie += NGHOST;
66  if (out_js != out_je) {
67  out_js -= NGHOST;
68  out_je += NGHOST;
69  }
70  if (out_ks != out_ke) {
71  out_ks -= NGHOST;
72  out_ke += NGHOST;
73  }
74  }
75 
76  // set ptrs to data in OutputData linked list, then slice/sum as needed
77  LoadOutputData(pmb);
78  if (!TransformOutputData(pmb)) {
79  ClearOutputData(); // required when LoadOutputData() is used.
80  continue;
81  } // skip if slice was out of range
82 
83  // create filename: "file_basename"+
84  // "."+"blockid"+"."+"file_id"+"."+XXXXX+".nc", where XXXXX = 5-digit
85  // file_number
86  std::string fname;
87  char number[6];
88  snprintf(number, sizeof(number), "%05d", output_params.file_number);
89  char blockid[12];
90  snprintf(blockid, sizeof(blockid), "block%d", pmb->gid);
91 
92  fname.assign(output_params.file_basename);
93  fname.append(".");
94  fname.append(blockid);
95  fname.append(".");
96  fname.append(output_params.file_id);
97  fname.append(".");
98  fname.append(number);
99  fname.append(".nc");
100 
101  // 1. open file for output
102  std::stringstream msg;
103  int ifile;
104 
105  nc_create(fname.c_str(), NC_NETCDF4, &ifile);
106 
107  // 2. coordinate structure
108  int ncells1 = out_ie - out_is + 1;
109  int ncells2 = out_je - out_js + 1;
110  int ncells3 = out_ke - out_ks + 1;
111  int nfaces1 = ncells1;
112  if (ncells1 > 1) nfaces1++;
113  int nfaces2 = ncells2;
114  if (ncells2 > 1) nfaces2++;
115  int nfaces3 = ncells3;
116  if (ncells3 > 1) nfaces3++;
117 
118  int nrays = prad->GetNumOutgoingRays();
119 
120  // 2. define coordinate
121  int idt, idx1, idx2, idx3, idx1f, idx2f, idx3f, iray;
122  // time
123  nc_def_dim(ifile, "time", NC_UNLIMITED, &idt);
124 
125  nc_def_dim(ifile, "x1", ncells1, &idx1);
126  if (ncells1 > 1) nc_def_dim(ifile, "x1f", nfaces1, &idx1f);
127 
128  nc_def_dim(ifile, "x2", ncells2, &idx2);
129  if (ncells2 > 1) nc_def_dim(ifile, "x2f", nfaces2, &idx2f);
130 
131  nc_def_dim(ifile, "x3", ncells3, &idx3);
132  if (ncells3 > 1) nc_def_dim(ifile, "x3f", nfaces3, &idx3f);
133 
134  if (nrays > 0) nc_def_dim(ifile, "ray", nrays, &iray);
135 
136  // 3. define variables
137  int ivt, ivx1, ivx2, ivx3, ivx1f, ivx2f, ivx3f, imu, iphi;
138  int loc[4] = {(int)pmb->loc.lx1, (int)pmb->loc.lx2, (int)pmb->loc.lx3,
139  pmb->loc.level};
140  int pos[4];
141 
142  nc_def_var(ifile, "time", NC_FLOAT, 1, &idt, &ivt);
143  nc_put_att_text(ifile, ivt, "axis", 1, "T");
144  nc_put_att_text(ifile, ivt, "units", 1, "s");
145  nc_put_att_text(ifile, ivt, "long_name", 4, "time");
146 
147  nc_def_var(ifile, "x1", NC_FLOAT, 1, &idx1, &ivx1);
148  nc_put_att_text(ifile, ivx1, "axis", 1, "Z");
149  nc_put_att_text(ifile, ivx1, "units", 1, "m");
150  nc_put_att_text(ifile, ivx1, "long_name", 27,
151  "Z-coordinate at cell center");
152 
153  pos[0] = 1;
154  pos[1] = pmb->pmy_mesh->mesh_size.nx1;
155  pos[2] = ncells1 * loc[0] + 1;
156  pos[3] = ncells1 * (loc[0] + 1);
157  nc_put_att_int(ifile, ivx1, "domain_decomposition", NC_INT, 4, pos);
158  if (ncells1 > 1) {
159  nc_def_var(ifile, "x1f", NC_FLOAT, 1, &idx1f, &ivx1f);
160  nc_put_att_text(ifile, ivx1f, "units", 1, "m");
161  nc_put_att_text(ifile, ivx1f, "long_name", 25,
162  "Z-coordinate at cell face");
163  pos[0]--;
164  pos[2]--;
165  nc_put_att_int(ifile, ivx1f, "domain_decomposition", NC_INT, 4, pos);
166  }
167 
168  nc_def_var(ifile, "x2", NC_FLOAT, 1, &idx2, &ivx2);
169  nc_put_att_text(ifile, ivx2, "axis", 1, "Y");
170  nc_put_att_text(ifile, ivx2, "units", 1, "m");
171  nc_put_att_text(ifile, ivx2, "long_name", 27,
172  "Y-coordinate at cell center");
173 
174  pos[0] = 1;
175  pos[1] = pmb->pmy_mesh->mesh_size.nx2;
176  pos[2] = ncells2 * loc[1] + 1;
177  pos[3] = ncells2 * (loc[1] + 1);
178  nc_put_att_int(ifile, ivx2, "domain_decomposition", NC_INT, 4, pos);
179  if (ncells2 > 1) {
180  nc_def_var(ifile, "x2f", NC_FLOAT, 1, &idx2f, &ivx2f);
181  nc_put_att_text(ifile, ivx2f, "units", 1, "m");
182  nc_put_att_text(ifile, ivx2f, "long_name", 25,
183  "Y-coordinate at cell face");
184  pos[0]--;
185  pos[2]--;
186  nc_put_att_int(ifile, ivx2f, "domain_decomposition", NC_INT, 4, pos);
187  }
188 
189  nc_def_var(ifile, "x3", NC_FLOAT, 1, &idx3, &ivx3);
190  nc_put_att_text(ifile, ivx3, "axis", 1, "X");
191  nc_put_att_text(ifile, ivx3, "units", 1, "m");
192  nc_put_att_text(ifile, ivx3, "long_name", 27,
193  "X-coordinate at cell center");
194 
195  pos[0] = 1;
196  pos[1] = pmb->pmy_mesh->mesh_size.nx3;
197  pos[2] = ncells3 * loc[2] + 1;
198  pos[3] = ncells3 * (loc[2] + 1);
199  nc_put_att_int(ifile, ivx3, "domain_decomposition", NC_INT, 4, pos);
200  if (ncells3 > 1) {
201  nc_def_var(ifile, "x3f", NC_FLOAT, 1, &idx3f, &ivx3f);
202  nc_put_att_text(ifile, ivx3f, "units", 1, "m");
203  nc_put_att_text(ifile, ivx3f, "long_name", 25,
204  "X-coordinate at cell face");
205  pos[0]--;
206  pos[2]--;
207  nc_put_att_int(ifile, ivx3f, "domain_decomposition", NC_INT, 4, pos);
208  }
209 
210  if (nrays > 0) {
211  nc_def_var(ifile, "mu_out", NC_FLOAT, 1, &iray, &imu);
212  nc_put_att_text(ifile, imu, "units", 1, "1");
213  nc_put_att_text(ifile, imu, "long_name", 18, "cosine polar angle");
214  nc_def_var(ifile, "phi_out", NC_FLOAT, 1, &iray, &iphi);
215  nc_put_att_text(ifile, iphi, "units", 3, "rad");
216  nc_put_att_text(ifile, iphi, "long_name", 15, "azimuthal angle");
217  }
218 
219  nc_put_att_int(ifile, NC_GLOBAL, "NumFilesInSet", NC_INT, 1, &pm->nbtotal);
220  if ((output_params.variable.compare("diag") == 0) ||
221  (output_params.variable.compare("flux") == 0)) {
222  float dt = (float)pm->dt / 2;
223  nc_put_att_float(ifile, NC_GLOBAL, "Time_shift_for_DiagFlux_variables",
224  NC_FLOAT, 1, &dt);
225  }
226 
227  OutputData *pdata = pfirst_data_;
228 
229  // count total variables (vector variables are expanded into flat scalars)
230  int total_vars = 0;
231  while (pdata != nullptr) {
232  std::string grid = pmeta->GetGridType(pdata->name);
233  int nvar = get_num_variables(grid, pdata->data);
234  total_vars += nvar;
235  pdata = pdata->pnext;
236  }
237 
238  int iaxis[4] = {idt, idx1, idx2, idx3};
239  int iaxis1[4] = {idt, idx1f, idx2, idx3};
240  int iaxis2[4] = {idt, idx1, idx2f, idx3};
241  int iaxis3[4] = {idt, idx1, idx2, idx3f};
242  int iaxisr[4] = {idt, iray, idx2, idx3};
243  int *var_ids = new int[total_vars];
244  int *ivar = var_ids;
245 
246  pdata = pfirst_data_;
247  while (pdata != nullptr) {
248  std::string name, attr;
249  std::vector<std::string> varnames;
250  std::string grid = pmeta->GetGridType(pdata->name);
251  int nvar = get_num_variables(grid, pdata->data);
252 
253  for (int n = 0; n < nvar; ++n) {
254  size_t pos = pdata->name.find('?');
255  if (nvar == 1) { // SCALARS
256  if (pos < pdata->name.length()) { // find '?'
257  varnames.push_back(pdata->name.substr(0, pos) +
258  pdata->name.substr(pos + 1));
259  } else {
260  varnames.push_back(pdata->name);
261  }
262  } else { // VECTORS
263  char c[16];
264  snprintf(c, sizeof(c), "%d", n + 1);
265  if (pos < pdata->name.length()) { // find '?'
266  varnames.push_back(pdata->name.substr(0, pos) + c +
267  pdata->name.substr(pos + 1));
268  } else {
269  varnames.push_back(pdata->name + c);
270  }
271  }
272  }
273 
274  for (int n = 0; n < nvar; ++n) {
275  name = varnames[n];
276  if (grid == "RCC") // radiation rays
277  nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxisr, ivar);
278  else if (grid == "CCF")
279  nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis1, ivar);
280  else if ((grid == "CFC") && (ncells2 > 1))
281  nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis2, ivar);
282  else if ((grid == "FCC") && (ncells3 > 1))
283  nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis3, ivar);
284  else if (grid == "--C")
285  nc_def_var(ifile, name.c_str(), NC_FLOAT, 2, iaxis, ivar);
286  else if (grid == "--F")
287  nc_def_var(ifile, name.c_str(), NC_FLOAT, 2, iaxis1, ivar);
288  else if (grid == "---")
289  nc_def_var(ifile, name.c_str(), NC_FLOAT, 1, iaxis, ivar);
290  else
291  nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis, ivar);
292 
293  // set units
294  attr = pmeta->GetUnits(pdata->name);
295  if (attr != "") {
296  nc_put_att_text(ifile, *ivar, "units", attr.length(), attr.c_str());
297  }
298 
299  // set long_name
300  attr = pmeta->GetLongName(pdata->name);
301  if (attr != "") {
302  nc_put_att_text(ifile, *ivar, "long_name", attr.length(),
303  attr.c_str());
304  }
305 
306  ivar++;
307  }
308  pdata = pdata->pnext;
309  }
310 
311  nc_enddef(ifile);
312 
313  // 4. write variables
314  float *data = new float[nfaces1 * nfaces2 * nfaces3];
315  size_t start[4] = {0, 0, 0, 0};
316  size_t count[4] = {1, (size_t)ncells1, (size_t)ncells2, (size_t)ncells3};
317  size_t count1[4] = {1, (size_t)nfaces1, (size_t)ncells2, (size_t)ncells3};
318  size_t count2[4] = {1, (size_t)ncells1, (size_t)nfaces2, (size_t)ncells3};
319  size_t count3[4] = {1, (size_t)ncells1, (size_t)ncells2, (size_t)nfaces3};
320  size_t countr[4] = {1, (size_t)nrays, (size_t)ncells2, (size_t)ncells3};
321 
322  float time = (float)pm->time;
323  nc_put_vara_float(ifile, ivt, start, count, &time);
324 
325  for (int i = out_is; i <= out_ie; ++i)
326  data[i - out_is] = (float)(pmb->pcoord->x1v(i));
327  nc_put_var_float(ifile, ivx1, data);
328 
329  if (ncells1 > 1) {
330  for (int i = out_is; i <= out_ie + 1; ++i)
331  data[i - out_is] = (float)(pmb->pcoord->x1f(i));
332  nc_put_var_float(ifile, ivx1f, data);
333  }
334 
335  for (int j = out_js; j <= out_je; ++j)
336  data[j - out_js] = (float)(pmb->pcoord->x2v(j));
337  nc_put_var_float(ifile, ivx2, data);
338 
339  if (ncells2 > 1) {
340  for (int j = out_js; j <= out_je + 1; ++j)
341  data[j - out_js] = (float)(pmb->pcoord->x2f(j));
342  nc_put_var_float(ifile, ivx2f, data);
343  }
344 
345  for (int k = out_ks; k <= out_ke; ++k)
346  data[k - out_ks] = (float)(pmb->pcoord->x3v(k));
347  nc_put_var_float(ifile, ivx3, data);
348 
349  if (ncells3 > 1) {
350  for (int k = out_ks; k <= out_ke + 1; ++k)
351  data[k - out_ks] = (float)(pmb->pcoord->x1f(k));
352  nc_put_var_float(ifile, ivx3f, data);
353  }
354 
355  if (nrays > 0) {
356  int m = 0;
357  for (int b = 0; b < prad->GetNumBands(); ++b) {
358  auto p = prad->GetBand(b);
359  for (int n = 0; n < p->GetNumOutgoingRays(); ++n)
360  data[m++] = (float)(p->GetCosinePolarAngle(n));
361  }
362  nc_put_var_float(ifile, imu, data);
363 
364  m = 0;
365  for (int b = 0; b < prad->GetNumBands(); ++b) {
366  auto p = prad->GetBand(b);
367  for (int n = 0; n < p->GetNumOutgoingRays(); ++n)
368  data[m++] = (float)(p->GetAzimuthalAngle(n));
369  }
370  nc_put_var_float(ifile, iphi, data);
371  }
372 
373  ivar = var_ids;
374  pdata = pfirst_data_;
375  while (pdata != nullptr) {
376  std::string grid = pmeta->GetGridType(pdata->name);
377  int nvar = get_num_variables(grid, pdata->data);
378 
379  if (grid == "RCC") { // radiation rays
380  for (int n = 0; n < nvar; n++) {
381  float *it = data;
382  for (int m = 0; m < nrays; ++m)
383  for (int j = out_js; j <= out_je; ++j)
384  for (int k = out_ks; k <= out_ke; ++k)
385  *it++ = (float)pdata->data(n, m, k, j);
386  nc_put_vara_float(ifile, *ivar++, start, countr, data);
387  }
388  } else if (grid == "CCF") {
389  for (int n = 0; n < nvar; n++) {
390  float *it = data;
391  for (int i = out_is; i <= out_ie + 1; ++i)
392  for (int j = out_js; j <= out_je; ++j)
393  for (int k = out_ks; k <= out_ke; ++k)
394  *it++ = (float)pdata->data(n, k, j, i);
395  nc_put_vara_float(ifile, *ivar++, start, count1, data);
396  }
397  } else if ((grid == "CFC") && (ncells2 > 1)) {
398  for (int n = 0; n < nvar; n++) {
399  float *it = data;
400  for (int i = out_is; i <= out_ie; ++i)
401  for (int j = out_js; j <= out_je + 1; ++j)
402  for (int k = out_ks; k <= out_ke; ++k)
403  *it++ = (float)pdata->data(n, k, j, i);
404  nc_put_vara_float(ifile, *ivar++, start, count2, data);
405  }
406  } else if ((grid == "FCC") && (ncells3 > 1)) {
407  for (int n = 0; n < nvar; n++) {
408  float *it = data;
409  for (int i = out_is; i <= out_ie; ++i)
410  for (int j = out_js; j <= out_je; ++j)
411  for (int k = out_ks; k <= out_ke + 1; ++k)
412  *it++ = (float)pdata->data(n, k, j, i);
413  nc_put_vara_float(ifile, *ivar++, start, count3, data);
414  }
415  } else if (grid == "--C") {
416  for (int n = 0; n < nvar; n++) {
417  float *it = data;
418  for (int i = out_is; i <= out_ie; ++i)
419  *it++ = (float)pdata->data(n, i);
420  nc_put_vara_float(ifile, *ivar++, start, count, data);
421  }
422  } else if (grid == "--F") {
423  for (int n = 0; n < nvar; n++) {
424  float *it = data;
425  for (int i = out_is; i <= out_ie + 1; ++i)
426  *it++ = (float)pdata->data(n, i);
427  nc_put_vara_float(ifile, *ivar++, start, count1, data);
428  }
429  } else if (grid == "---") {
430  for (int n = 0; n < nvar; n++) {
431  float *it = data;
432  *it++ = (float)pdata->data(n);
433  nc_put_vara_float(ifile, *ivar++, start, count, data);
434  }
435  } else {
436  for (int n = 0; n < nvar; n++) {
437  float *it = data;
438  for (int i = out_is; i <= out_ie; ++i)
439  for (int j = out_js; j <= out_je; ++j)
440  for (int k = out_ks; k <= out_ke; ++k)
441  *it++ = (float)pdata->data(n, k, j, i);
442  nc_put_vara_float(ifile, *ivar++, start, count, data);
443  }
444  }
445 
446  // doesn't work
447  // nc_put_att_text(ifile, *(ivar-1), "output",
448  // output_params.variable.length(), output_params.variable.c_str());
449  pdata = pdata->pnext;
450  }
451 
452  // 5. close nc file
453  nc_close(ifile);
454 
455  ClearOutputData(); // required when LoadOutputData() is used.
456  delete[] data;
457  delete[] var_ids;
458  } // end loop over MeshBlocks
459 
460  // increment counters
461  output_params.file_number++;
462  output_params.next_time += output_params.dt;
463  pin->SetInteger(output_params.block_name, "file_number",
464  output_params.file_number);
465  pin->SetReal(output_params.block_name, "next_time", output_params.next_time);
466 
467  return;
468 }
469 
470 #endif // NETCDFOUTPUT
MetadataTable::GetInstance
static MetadataTable const * GetInstance()
Definition: output_utils.cpp:57
fit_ammonia.data
data
Definition: fit_ammonia.py:5
fit_ammonia.b
b
Definition: fit_ammonia.py:8
fit_ammonia.p
p
Definition: fit_ammonia.py:8
make_plots.j
j
Definition: make_plots.py:26
make_plots.time
time
Definition: make_plots.py:9
get_num_variables
int get_num_variables(std::string grid, AthenaArray< Real > const &data)
Definition: output_utils.cpp:12
output_utils.hpp
radiation_band.hpp