PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
NC_Serial.cc
Go to the documentation of this file.
1 // Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2019, 2020, 2023 PISM Authors
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 
19 #include <mpi.h>
20 
21 #include "pism/util/io/NC_Serial.hh"
22 
23 // The following is a stupid kludge necessary to make NetCDF 4.x work in
24 // serial mode in an MPI program:
25 #ifndef MPI_INCLUDED
26 #define MPI_INCLUDED 1
27 #endif
28 #include <netcdf.h>
29 
30 #include <cstdio> // stderr, fprintf
31 
32 #include "pism/util/pism_utilities.hh" // join
33 #include "pism/util/error_handling.hh"
34 
35 #include "pism/util/io/pism_type_conversion.hh" // This has to be included *after* netcdf.h.
36 
37 namespace pism {
38 namespace io {
39 
40 //! \brief Prints an error message; for debugging.
41 static void check(const ErrorLocation &where, int return_code) {
42  if (return_code != NC_NOERR) {
43  throw RuntimeError(where, nc_strerror(return_code));
44  }
45 }
46 
47 //! call MPI_Abort() if a NetCDF call failed
48 static void check_and_abort(MPI_Comm com, const ErrorLocation &where, int return_code) {
49  if (return_code != NC_NOERR) {
50  fprintf(stderr, "%s:%d: %s\n", where.filename, where.line_number, nc_strerror(return_code));
51  MPI_Abort(com, -1);
52  }
53 }
54 
56  : NCFile(c), m_rank(0) {
57  MPI_Comm_rank(m_com, &m_rank);
58 }
59 
61  if (m_file_id >= 0) {
62  if (m_rank == 0) {
63  nc_close(m_file_id);
64  fprintf(stderr, "NC_Serial::~NC_Serial: NetCDF file %s is still open\n",
65  m_filename.c_str());
66  }
67  m_file_id = -1;
68  }
69 }
70 
72  (void) level;
73  // NetCDF-3 does not support compression.
74 }
75 
76 // open/create/close
77 void NC_Serial::open_impl(const std::string &fname, io::Mode mode) {
78  int stat = NC_NOERR;
79 
80  int open_mode = mode == io::PISM_READONLY ? NC_NOWRITE : NC_WRITE;
81 
82  if (m_rank == 0) {
83  stat = nc_open(fname.c_str(), open_mode, &m_file_id);
84  }
85 
86  MPI_Barrier(m_com);
87  MPI_Bcast(&m_file_id, 1, MPI_INT, 0, m_com);
88  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
89 
91 }
92 
93 //! \brief Create a NetCDF file.
94 void NC_Serial::create_impl(const std::string &fname) {
95  int stat = NC_NOERR;
96 
97  if (m_rank == 0) {
98  stat = nc_create(fname.c_str(), NC_CLOBBER | NC_64BIT_OFFSET, &m_file_id);
99  }
100 
101  MPI_Barrier(m_com);
102  MPI_Bcast(&m_file_id, 1, MPI_INT, 0, m_com);
103  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
104 
105  check(PISM_ERROR_LOCATION, stat);
106 }
107 
108 //! \brief Close a NetCDF file.
110  int stat = NC_NOERR;
111 
112  if (m_rank == 0) {
113  stat = nc_close(m_file_id);
114  }
115 
116  m_file_id = -1;
117 
118  MPI_Barrier(m_com);
119  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
120 
121  check(PISM_ERROR_LOCATION, stat);
122 }
123 
124 
125 void NC_Serial::sync_impl() const {
126  int stat = NC_NOERR;
127 
128  if (m_rank == 0) {
129  stat = nc_sync(m_file_id);
130  }
131 
132  MPI_Barrier(m_com);
133  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
134 
135  check(PISM_ERROR_LOCATION, stat);
136 }
137 
138 
139 //! \brief Exit define mode.
141  int stat = NC_NOERR;
142 
143  int header_size = 200 * 1024;
144 
145  if (m_rank == 0) {
146  stat = nc__enddef(m_file_id, header_size, 4, 0, 4);
147  }
148 
149  MPI_Barrier(m_com);
150  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
151 
152  check(PISM_ERROR_LOCATION, stat);
153 }
154 
155 //! \brief Enter define mode.
156 void NC_Serial::redef_impl() const {
157  int stat = NC_NOERR;
158 
159  if (m_rank == 0) {
160  stat = nc_redef(m_file_id);
161  }
162 
163  MPI_Barrier(m_com);
164  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
165 
166  check(PISM_ERROR_LOCATION, stat);
167 }
168 
169 
170 //! \brief Define a dimension.
171 void NC_Serial::def_dim_impl(const std::string &name, size_t length) const {
172  int stat = NC_NOERR;
173 
174  if (m_rank == 0) {
175  int dimid;
176  stat = nc_def_dim(m_file_id, name.c_str(), length, &dimid);
177  }
178 
179  MPI_Barrier(m_com);
180  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
181 
182  check(PISM_ERROR_LOCATION, stat);
183 }
184 
185 void NC_Serial::inq_dimid_impl(const std::string &dimension_name, bool &exists) const {
186  int stat, flag = -1;
187 
188  if (m_rank == 0) {
189  stat = nc_inq_dimid(m_file_id, dimension_name.c_str(), &flag);
190 
191  flag = (stat == NC_NOERR) ? 1 : 0;
192  }
193  MPI_Barrier(m_com);
194  MPI_Bcast(&flag, 1, MPI_INT, 0, m_com);
195 
196  exists = (flag == 1);
197 }
198 
199 
200 //! \brief Get a dimension length.
201 void NC_Serial::inq_dimlen_impl(const std::string &dimension_name, unsigned int &result) const {
202  int stat = NC_NOERR;
203 
204  if (m_rank == 0) {
205  int dimid;
206  size_t length;
207 
208  stat = nc_inq_dimid(m_file_id, dimension_name.c_str(), &dimid);
209 
210  if (stat == NC_NOERR) {
211  stat = nc_inq_dimlen(m_file_id, dimid, &length);
212  result = static_cast<unsigned int>(length);
213  }
214  }
215 
216  MPI_Barrier(m_com);
217  MPI_Bcast(&result, 1, MPI_UNSIGNED, 0, m_com);
218  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
219 
220  check(PISM_ERROR_LOCATION, stat);
221 }
222 
223 //! \brief Get an unlimited dimension.
224 void NC_Serial::inq_unlimdim_impl(std::string &result) const {
225  int stat = NC_NOERR;
226  std::vector<char> dimname(NC_MAX_NAME + 1, 0);
227 
228  if (m_rank == 0) {
229  int dimid;
230  stat = nc_inq_unlimdim(m_file_id, &dimid);
231 
232  // nc_inq_unlimdim() sets dimid to -1 if there is no unlimited dimension
233  if (stat == NC_NOERR and dimid != -1) {
234  stat = nc_inq_dimname(m_file_id, dimid, dimname.data());
235  } else {
236  // leave dimname empty
237  stat = NC_NOERR;
238  }
239  }
240 
241  MPI_Barrier(m_com);
242 
243  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
244  MPI_Bcast(dimname.data(), NC_MAX_NAME, MPI_CHAR, 0, m_com);
245 
246  check(PISM_ERROR_LOCATION, stat);
247 
248  result = dimname.data();
249 }
250 
251 //! \brief Define a variable.
252 void NC_Serial::def_var_impl(const std::string &name, io::Type nctype,
253  const std::vector<std::string> &dims) const {
254  int stat = NC_NOERR;
255 
256  if (m_rank == 0) {
257  std::vector<int> dimids;
258  int varid;
259 
260  for (auto d : dims) {
261  int dimid;
262  stat = nc_inq_dimid(m_file_id, d.c_str(), &dimid);
263  if (stat != NC_NOERR) {
264  break;
265  }
266  dimids.push_back(dimid);
267  }
268 
269  if (stat == NC_NOERR) {
270  stat = nc_def_var(m_file_id, name.c_str(), pism_type_to_nc_type(nctype),
271  static_cast<int>(dims.size()), dimids.data(), &varid);
272  }
273  }
274 
275  MPI_Barrier(m_com);
276  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
277 
278  check(PISM_ERROR_LOCATION, stat);
279 }
280 
281 void NC_Serial::get_varm_double_impl(const std::string &variable_name,
282  const std::vector<unsigned int> &start,
283  const std::vector<unsigned int> &count,
284  const std::vector<unsigned int> &imap, double *op) const {
285  return this->get_var_double(variable_name, start, count, imap, op, true);
286 }
287 
288 void NC_Serial::get_vara_double_impl(const std::string &variable_name,
289  const std::vector<unsigned int> &start,
290  const std::vector<unsigned int> &count, double *op) const {
291  std::vector<unsigned int> dummy;
292  return this->get_var_double(variable_name, start, count, dummy, op, false);
293 }
294 
295 //! \brief Get variable data.
296 void NC_Serial::get_var_double(const std::string &variable_name,
297  const std::vector<unsigned int> &start_input,
298  const std::vector<unsigned int> &count_input,
299  const std::vector<unsigned int> &imap_input, double *ip,
300  bool transposed) const {
301  std::vector<unsigned int> start = start_input;
302  std::vector<unsigned int> count = count_input;
303  std::vector<unsigned int> imap = imap_input;
304  const int start_tag = 1, count_tag = 2, data_tag = 3, imap_tag = 4, chunk_size_tag = 5;
305  int stat = NC_NOERR, com_size, ndims = static_cast<int>(start.size());
306  MPI_Status mpi_stat;
307  unsigned int local_chunk_size = 1, processor_0_chunk_size = 0;
308 
309  if (not transposed) {
310  imap.resize(ndims);
311  }
312 
313  // get the size of the communicator
314  MPI_Comm_size(m_com, &com_size);
315 
316  // compute the size of a local chunk
317  for (int k = 0; k < ndims; ++k) {
318  local_chunk_size *= count[k];
319  }
320 
321  // compute the maximum and send it to processor 0; this is the size of the
322  // buffer processor 0 will need
323  MPI_Reduce(&local_chunk_size, &processor_0_chunk_size, 1, MPI_UNSIGNED, MPI_MAX, 0, m_com);
324 
325  // now we need to send start, count and imap data to processor 0 and receive data
326  if (m_rank == 0) {
327  std::vector<double> processor_0_buffer;
328  // Note: this could be optimized: if processor_0_chunk_size <=
329  // max(local_chunk_size) we can avoid allocating this buffer. The inner for
330  // loop will have to be re-ordered, though.
331  processor_0_buffer.resize(processor_0_chunk_size);
332 
333  // MPI calls below require C datatypes (so that we don't have to worry
334  // about sizes of size_t and ptrdiff_t), so we make local copies of start,
335  // count, and imap to use in the nc_get_varm_double() call.
336  std::vector<size_t> nc_start(ndims), nc_count(ndims);
337  std::vector<ptrdiff_t> nc_imap(ndims), nc_stride(ndims);
338  int varid;
339 
340  stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
342 
343  for (int r = 0; r < com_size; ++r) {
344 
345  if (r != 0) {
346  // Note: start, count, imap, and local_chunk_size on processor zero are
347  // used *before* they get overwritten by these calls
348  MPI_Recv(start.data(), ndims, MPI_UNSIGNED, r, start_tag, m_com, &mpi_stat);
349  MPI_Recv(count.data(), ndims, MPI_UNSIGNED, r, count_tag, m_com, &mpi_stat);
350  MPI_Recv(imap.data(), ndims, MPI_UNSIGNED, r, imap_tag, m_com, &mpi_stat);
351  MPI_Recv(&local_chunk_size, 1, MPI_UNSIGNED, r, chunk_size_tag, m_com, &mpi_stat);
352  }
353 
354  // This for loop uses start, count and imap passed in as arguments when r
355  // == 0. For r > 0 they are overwritten by MPI_Recv calls above.
356  for (int k = 0; k < ndims; ++k) {
357  nc_start[k] = start[k];
358  nc_count[k] = count[k];
359  nc_imap[k] = imap[k];
360  nc_stride[k] = 1; // fill with ones; this way it works even with
361  // NetCDF versions with a bug affecting the
362  // stride == NULL case.
363  }
364 
365  if (transposed) {
366  stat = nc_get_varm_double(m_file_id, varid, nc_start.data(), nc_count.data(),
367  nc_stride.data(), nc_imap.data(), processor_0_buffer.data());
368  } else {
369  stat = nc_get_vara_double(m_file_id, varid, nc_start.data(), nc_count.data(),
370  processor_0_buffer.data());
371  }
372  check(PISM_ERROR_LOCATION, stat);
373 
374  if (r != 0) {
375  MPI_Send(processor_0_buffer.data(), (int)local_chunk_size, MPI_DOUBLE, r, data_tag, m_com);
376  } else {
377  for (unsigned int k = 0; k < local_chunk_size; ++k) {
378  ip[k] = processor_0_buffer[k];
379  }
380  }
381  } // end of the for loop
382 
383  } else {
384  MPI_Send(start.data(), ndims, MPI_UNSIGNED, 0, start_tag, m_com);
385  MPI_Send(count.data(), ndims, MPI_UNSIGNED, 0, count_tag, m_com);
386  MPI_Send(imap.data(), ndims, MPI_UNSIGNED, 0, imap_tag, m_com);
387  MPI_Send(&local_chunk_size, 1, MPI_UNSIGNED, 0, chunk_size_tag, m_com);
388 
389  MPI_Recv(ip, (int)local_chunk_size, MPI_DOUBLE, 0, data_tag, m_com, &mpi_stat);
390  }
391 }
392 
393 void NC_Serial::put_vara_double_impl(const std::string &variable_name,
394  const std::vector<unsigned int> &start_input,
395  const std::vector<unsigned int> &count_input,
396  const double *op) const {
397  // make copies of start and count so that we can use them in MPI_Recv() calls below
398  std::vector<unsigned int> start = start_input;
399  std::vector<unsigned int> count = count_input;
400  const int start_tag = 1, count_tag = 2, data_tag = 3, chunk_size_tag = 4;
401  int stat = NC_NOERR, com_size = 0, ndims = static_cast<int>(start.size());
402  MPI_Status mpi_stat;
403  unsigned int local_chunk_size = 1, processor_0_chunk_size = 0;
404 
405  // get the size of the communicator
406  MPI_Comm_size(m_com, &com_size);
407 
408  // compute the size of a local chunk
409  for (int k = 0; k < ndims; ++k) {
410  local_chunk_size *= count[k];
411  }
412 
413  // compute the maximum and send it to processor 0; this is the size of the
414  // buffer processor 0 will need
415  MPI_Reduce(&local_chunk_size, &processor_0_chunk_size, 1, MPI_UNSIGNED, MPI_MAX, 0, m_com);
416 
417  // now we need to send start and count data to processor 0 and receive data
418  if (m_rank == 0) {
419  std::vector<double> processor_0_buffer;
420  processor_0_buffer.resize(processor_0_chunk_size);
421 
422  // MPI calls below require C datatypes (so that we don't have to worry about sizes of
423  // size_t and ptrdiff_t), so we make local copies of start and count to use in the
424  // nc_get_vara_double() call.
425  std::vector<size_t> nc_start(ndims), nc_count(ndims);
426  std::vector<ptrdiff_t> nc_stride(ndims);
427  int varid;
428 
429  stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
431 
432  for (int r = 0; r < com_size; ++r) {
433 
434  if (r != 0) {
435  // Note: start, count, and local_chunk_size on processor zero are used *before*
436  // they get overwritten by these calls
437  MPI_Recv(start.data(), ndims, MPI_UNSIGNED, r, start_tag, m_com, &mpi_stat);
438  MPI_Recv(count.data(), ndims, MPI_UNSIGNED, r, count_tag, m_com, &mpi_stat);
439  MPI_Recv(&local_chunk_size, 1, MPI_UNSIGNED, r, chunk_size_tag, m_com, &mpi_stat);
440 
441  MPI_Recv(processor_0_buffer.data(), local_chunk_size, MPI_DOUBLE, r, data_tag, m_com,
442  &mpi_stat);
443  } else {
444  for (unsigned int k = 0; k < local_chunk_size; ++k) {
445  processor_0_buffer[k] = op[k];
446  }
447  }
448 
449  // This for loop uses start and count passed in as arguments when r == 0. For r > 0
450  // they are overwritten by MPI_Recv calls above.
451  for (int k = 0; k < ndims; ++k) {
452  nc_start[k] = start[k];
453  nc_count[k] = count[k];
454  nc_stride[k] = 1; // fill with ones; this way it works even with
455  // NetCDF versions with a bug affecting the
456  // stride == NULL case.
457  }
458 
459  stat = nc_put_vara_double(m_file_id, varid, nc_start.data(), nc_count.data(),
460  processor_0_buffer.data());
461  check(PISM_ERROR_LOCATION, stat);
462  } // end of the for loop
463  } else {
464  MPI_Send(start.data(), ndims, MPI_UNSIGNED, 0, start_tag, m_com);
465  MPI_Send(count.data(), ndims, MPI_UNSIGNED, 0, count_tag, m_com);
466  MPI_Send(&local_chunk_size, 1, MPI_UNSIGNED, 0, chunk_size_tag, m_com);
467 
468  MPI_Send(const_cast<double *>(op), (int)local_chunk_size, MPI_DOUBLE, 0, data_tag, m_com);
469  }
470 }
471 
472 //! \brief Get the number of variables.
473 void NC_Serial::inq_nvars_impl(int &result) const {
474  int stat = NC_NOERR;
475 
476  if (m_rank == 0) {
477  stat = nc_inq_nvars(m_file_id, &result);
478  }
479  MPI_Barrier(m_com);
480 
481  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
482  check(PISM_ERROR_LOCATION, stat);
483 
484  MPI_Bcast(&result, 1, MPI_INT, 0, m_com);
485 }
486 
487 //! \brief Get dimensions a variable depends on.
488 void NC_Serial::inq_vardimid_impl(const std::string &variable_name,
489  std::vector<std::string> &result) const {
490  int stat, ndims, varid = -1;
491  std::vector<int> dimids;
492 
493  if (m_rank == 0) {
494  stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
495 
496  if (stat == NC_NOERR) {
497  stat = nc_inq_varndims(m_file_id, varid, &ndims);
498  }
499  }
500 
501  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
502  check(PISM_ERROR_LOCATION, stat);
503 
504  MPI_Bcast(&ndims, 1, MPI_INT, 0, m_com);
505 
506  if (ndims == 0) {
507  result.clear();
508  return;
509  }
510 
511  result.resize(ndims);
512  dimids.resize(ndims);
513 
514  if (m_rank == 0) {
515  stat = nc_inq_vardimid(m_file_id, varid, dimids.data());
516  }
517 
518  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
519  check(PISM_ERROR_LOCATION, stat);
520 
521  MPI_Barrier(m_com);
522 
523  for (int k = 0; k < ndims; ++k) {
524  std::vector<char> name(NC_MAX_NAME + 1, 0);
525 
526  if (m_rank == 0) {
527  stat = nc_inq_dimname(m_file_id, dimids[k], name.data());
528  }
529 
530  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
531  check(PISM_ERROR_LOCATION, stat);
532 
533  MPI_Barrier(m_com);
534  MPI_Bcast(name.data(), name.size(), MPI_CHAR, 0, m_com);
535 
536  result[k] = name.data();
537  }
538 }
539 
540 //! \brief Get the number of attributes of a variable.
541 /*!
542  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
543  */
544 void NC_Serial::inq_varnatts_impl(const std::string &variable_name, int &result) const {
545  int stat = NC_NOERR;
546 
547  if (m_rank == 0) {
548  int varid = get_varid(variable_name);
549 
550  if (varid >= NC_GLOBAL) {
551  stat = nc_inq_varnatts(m_file_id, varid, &result);
552  } else {
553  stat = varid; // LCOV_EXCL_LINE
554  }
555  }
556  MPI_Barrier(m_com);
557 
558  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
559  check(PISM_ERROR_LOCATION, stat);
560 
561  MPI_Bcast(&result, 1, MPI_INT, 0, m_com);
562 }
563 
564 //! \brief Finds a variable and sets the "exists" flag.
565 void NC_Serial::inq_varid_impl(const std::string &variable_name, bool &exists) const {
566  int stat, flag = -1;
567 
568  if (m_rank == 0) {
569  stat = nc_inq_varid(m_file_id, variable_name.c_str(), &flag);
570 
571  flag = (stat == NC_NOERR) ? 1 : 0;
572  }
573  MPI_Barrier(m_com);
574  MPI_Bcast(&flag, 1, MPI_INT, 0, m_com);
575 
576  exists = (flag == 1);
577 }
578 
579 void NC_Serial::inq_varname_impl(unsigned int j, std::string &result) const {
580  int stat = NC_NOERR;
581  std::vector<char> varname(NC_MAX_NAME + 1, 0);
582 
583  if (m_rank == 0) {
584  stat = nc_inq_varname(m_file_id, j, varname.data());
585  }
586 
587  MPI_Barrier(m_com);
588 
589  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
590  MPI_Bcast(varname.data(), NC_MAX_NAME, MPI_CHAR, 0, m_com);
591 
592  check(PISM_ERROR_LOCATION, stat);
593 
594  result = varname.data();
595 }
596 
597 //! \brief Gets a double attribute.
598 /*!
599  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
600  */
601 void NC_Serial::get_att_double_impl(const std::string &variable_name, const std::string &att_name,
602  std::vector<double> &result) const {
603  int stat = NC_NOERR, len = 0;
604 
605  int varid = get_varid(variable_name);
606 
607  // Read and broadcast the attribute length:
608  if (m_rank == 0) {
609  size_t attlen = 0;
610 
611  if (varid >= NC_GLOBAL) {
612  stat = nc_inq_attlen(m_file_id, varid, att_name.c_str(), &attlen);
613  } else {
614  stat = varid; // LCOV_EXCL_LINE
615  }
616 
617  if (stat == NC_NOERR) {
618  len = static_cast<int>(attlen);
619  } else {
620  len = 0;
621  }
622  }
623  MPI_Bcast(&len, 1, MPI_INT, 0, m_com);
624 
625  if (len == 0) {
626  result.clear();
627  return;
628  }
629 
630  result.resize(len);
631 
632  // Now read data and broadcast stat to see if we succeeded:
633  if (m_rank == 0) {
634  stat = nc_get_att_double(m_file_id, varid, att_name.c_str(), result.data());
635  }
636  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
637 
638  check(PISM_ERROR_LOCATION, stat);
639 
640  // Broadcast data
641  MPI_Bcast(result.data(), len, MPI_DOUBLE, 0, m_com);
642 }
643 
644 // Get a text (character array) attribute on rank 0.
645 static int get_att_text(int ncid, int varid, const std::string &att_name, std::string &result) {
646  int stat = NC_NOERR;
647 
648  size_t attlen = 0;
649  stat = nc_inq_attlen(ncid, varid, att_name.c_str(), &attlen);
650  if (stat != NC_NOERR) {
651  result = "";
652  return 0;
653  }
654 
655  std::vector<char> buffer(attlen + 1, 0);
656  stat = nc_get_att_text(ncid, varid, att_name.c_str(), buffer.data());
657  if (stat == NC_NOERR) {
658  result = buffer.data();
659  } else {
660  result = "";
661  }
662 
663  return 0;
664 }
665 
666 // Get a string attribute on rank 0. In "string array" attributes array elements are concatenated
667 // using "," as the separator.
668 static int get_att_string(int ncid, int varid, const std::string &att_name, std::string &result) {
669  int stat = NC_NOERR;
670 
671  size_t attlen = 0;
672  stat = nc_inq_attlen(ncid, varid, att_name.c_str(), &attlen);
673  if (stat != NC_NOERR) {
674  result = "";
675  return 0;
676  }
677 
678  std::vector<char *> buffer(attlen + 1, 0);
679  stat = nc_get_att_string(ncid, varid, att_name.c_str(), buffer.data());
680  if (stat == NC_NOERR) {
681  std::vector<std::string> strings(attlen);
682  for (size_t k = 0; k < attlen; ++k) {
683  strings[k] = buffer[k];
684  }
685  result = join(strings, ",");
686  } else {
687  result = "";
688  }
689  stat = nc_free_string(attlen, buffer.data());
690 
691  return stat;
692 }
693 
694 
695 //! \brief Gets a text attribute.
696 /*!
697  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
698  */
699 void NC_Serial::get_att_text_impl(const std::string &variable_name, const std::string &att_name,
700  std::string &result) const {
701  int stat = NC_NOERR;
702 
703  // Read and broadcast the attribute length:
704  if (m_rank == 0) {
705 
706  int varid = get_varid(variable_name);
707 
708  if (varid >= NC_GLOBAL) {
709  nc_type nctype = NC_NAT;
710  stat = nc_inq_atttype(m_file_id, varid, att_name.c_str(), &nctype);
711 
712  if (stat == NC_NOERR) {
713  switch (nctype) {
714  case NC_CHAR:
715  stat = pism::io::get_att_text(m_file_id, varid, att_name, result);
716  break;
717  case NC_STRING:
718  stat = pism::io::get_att_string(m_file_id, varid, att_name, result);
719  break;
720  default:
721  result = "";
722  stat = NC_NOERR;
723  }
724  } else if (stat == NC_ENOTATT) {
725  result = "";
726  stat = NC_NOERR;
727  }
728  } else {
729  stat = varid; // LCOV_EXCL_LINE
730  }
731  }
732  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
733  check(PISM_ERROR_LOCATION, stat);
734 
735  unsigned int len = result.size();
736  MPI_Bcast(&len, 1, MPI_UNSIGNED, 0, m_com);
737 
738  result.resize(len);
739  MPI_Bcast(&result[0], (int)len, MPI_CHAR, 0, m_com);
740 }
741 
742 
743 //! \brief Writes a double attribute.
744 /*!
745  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
746  */
747 void NC_Serial::put_att_double_impl(const std::string &variable_name, const std::string &att_name,
748  io::Type nctype, const std::vector<double> &data) const {
749  int stat = NC_NOERR;
750 
751  if (m_rank == 0) {
752  int varid = get_varid(variable_name);
753 
754  if (varid >= NC_GLOBAL) {
755  stat = nc_put_att_double(m_file_id, varid, att_name.c_str(), pism_type_to_nc_type(nctype),
756  data.size(), data.data());
757  } else {
758  stat = varid; // LCOV_EXCL_LINE
759  }
760  }
761 
762  MPI_Barrier(m_com);
763  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
764 
765  check(PISM_ERROR_LOCATION, stat);
766 }
767 
768 
769 //! \brief Writes a text attribute.
770 /*!
771  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
772  */
773 void NC_Serial::put_att_text_impl(const std::string &variable_name, const std::string &att_name,
774  const std::string &value) const {
775  int stat = NC_NOERR;
776 
777  if (m_rank == 0) {
778  int varid = get_varid(variable_name);
779 
780  if (varid >= NC_GLOBAL) {
781  stat = nc_put_att_text(m_file_id, varid, att_name.c_str(), value.size(), value.c_str());
782  } else {
783  stat = varid; // LCOV_EXCL_LINE
784  }
785  }
786 
787  MPI_Barrier(m_com);
788  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
789 
790  check(PISM_ERROR_LOCATION, stat);
791 }
792 
793 //! \brief Gets the name of a numbered attribute.
794 /*!
795  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
796  */
797 void NC_Serial::inq_attname_impl(const std::string &variable_name, unsigned int n,
798  std::string &result) const {
799  int stat = NC_NOERR;
800  std::vector<char> name(NC_MAX_NAME + 1, 0);
801 
802  if (m_rank == 0) {
803  int varid = get_varid(variable_name);
804 
805  if (varid >= NC_GLOBAL) {
806  stat = nc_inq_attname(m_file_id, varid, n, name.data());
807  check(PISM_ERROR_LOCATION, stat);
808  } else {
809  stat = varid; // LCOV_EXCL_LINE
810  }
811  }
812  MPI_Barrier(m_com);
813  MPI_Bcast(name.data(), NC_MAX_NAME, MPI_CHAR, 0, m_com);
814  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
815 
816  check(PISM_ERROR_LOCATION, stat);
817 
818  result = name.data();
819 }
820 
821 //! \brief Gets the type of an attribute.
822 /*!
823  * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
824  */
825 void NC_Serial::inq_atttype_impl(const std::string &variable_name, const std::string &att_name,
826  io::Type &result) const {
827  int stat, tmp;
828 
829  if (m_rank == 0) {
830  int varid = get_varid(variable_name);
831 
832  if (varid >= NC_GLOBAL) {
833  // In NetCDF 3.6.x nc_type is an enum; in 4.x it is 'typedef int'.
834  nc_type nctype = NC_NAT;
835  stat = nc_inq_atttype(m_file_id, varid, att_name.c_str(), &nctype);
836  if (stat == NC_ENOTATT) {
837  tmp = NC_NAT;
838  stat = NC_NOERR;
839  } else {
840  tmp = static_cast<int>(nctype);
841  }
842  } else {
843  stat = varid; // LCOV_EXCL_LINE
844  }
845  }
846  MPI_Barrier(m_com);
847  MPI_Bcast(&tmp, 1, MPI_INT, 0, m_com);
848 
849  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
850  check(PISM_ERROR_LOCATION, stat);
851 
852  result = nc_type_to_pism_type(tmp);
853 }
854 
855 
856 //! \brief Sets the fill mode.
857 void NC_Serial::set_fill_impl(int fillmode, int &old_modep) const {
858  int stat = NC_NOERR;
859 
860  if (m_rank == 0) {
861  stat = nc_set_fill(m_file_id, fillmode, &old_modep);
862  }
863 
864  MPI_Barrier(m_com);
865  MPI_Bcast(&old_modep, 1, MPI_INT, 0, m_com);
866  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
867 
868  check(PISM_ERROR_LOCATION, stat);
869 }
870 
871 std::string NC_Serial::get_format() const {
872  int format;
873 
874  if (m_rank == 0) {
875  int stat = nc_inq_format(m_file_id, &format);
876  check(PISM_ERROR_LOCATION, stat);
877  }
878  MPI_Barrier(m_com);
879  MPI_Bcast(&format, 1, MPI_INT, 0, m_com);
880 
881  switch (format) {
882  case NC_FORMAT_CLASSIC:
883  case NC_FORMAT_64BIT_OFFSET:
884  return "netcdf3";
885  case NC_FORMAT_64BIT_DATA:
886  return "cdf5";
887  case NC_FORMAT_NETCDF4:
888  case NC_FORMAT_NETCDF4_CLASSIC:
889  default:
890  return "netcdf4";
891  }
892 }
893 
894 void NC_Serial::del_att_impl(const std::string &variable_name, const std::string &att_name) const {
895  int stat = NC_NOERR;
896 
897  if (m_rank == 0) {
898  int varid = get_varid(variable_name);
899 
900  if (varid >= NC_GLOBAL) {
901  stat = nc_del_att(m_file_id, varid, att_name.c_str());
902  }
903  }
904 
905  MPI_Barrier(m_com);
906  MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
907 
908  check(PISM_ERROR_LOCATION, stat);
909 }
910 
911 /*!
912  * return the varid corresponding to a variable.
913  *
914  * If the value returned is NC_GLOBAL or greater, it is a varid, otherwise it is an error
915  * code.
916  */
917 int NC_Serial::get_varid(const std::string &variable_name) const {
918  if (variable_name == "PISM_GLOBAL") {
919  return NC_GLOBAL;
920  }
921 
922  if (m_rank == 0) {
923  int varid = -2;
924  int stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
925 
926  if (stat == NC_NOERR) {
927  return varid;
928  }
929 
930  return stat;
931  }
932 
933  return -2; // this value will not be used
934 }
935 
936 } // end of namespace io
937 } // end of namespace pism
const char * filename
std::string m_filename
Definition: NCFile.hh:234
MPI_Comm m_com
Definition: NCFile.hh:232
The PISM wrapper for a subset of the NetCDF C API.
Definition: NCFile.hh:59
void get_att_text_impl(const std::string &variable_name, const std::string &att_name, std::string &result) const
Gets a text attribute.
Definition: NC_Serial.cc:699
virtual void create_impl(const std::string &filename)
Create a NetCDF file.
Definition: NC_Serial.cc:94
void def_dim_impl(const std::string &name, size_t length) const
Define a dimension.
Definition: NC_Serial.cc:171
void inq_varname_impl(unsigned int j, std::string &result) const
Definition: NC_Serial.cc:579
void put_vara_double_impl(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, const double *op) const
Definition: NC_Serial.cc:393
void set_fill_impl(int fillmode, int &old_modep) const
Sets the fill mode.
Definition: NC_Serial.cc:857
std::string get_format() const
Definition: NC_Serial.cc:871
void enddef_impl() const
Exit define mode.
Definition: NC_Serial.cc:140
void get_vara_double_impl(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, double *ip) const
Definition: NC_Serial.cc:288
void inq_vardimid_impl(const std::string &variable_name, std::vector< std::string > &result) const
Get dimensions a variable depends on.
Definition: NC_Serial.cc:488
void put_att_text_impl(const std::string &variable_name, const std::string &att_name, const std::string &value) const
Writes a text attribute.
Definition: NC_Serial.cc:773
void get_att_double_impl(const std::string &variable_name, const std::string &att_name, std::vector< double > &result) const
Gets a double attribute.
Definition: NC_Serial.cc:601
void get_varm_double_impl(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, const std::vector< unsigned int > &imap, double *ip) const
Definition: NC_Serial.cc:281
void put_att_double_impl(const std::string &variable_name, const std::string &att_name, io::Type xtype, const std::vector< double > &data) const
Writes a double attribute.
Definition: NC_Serial.cc:747
void open_impl(const std::string &filename, io::Mode mode)
Definition: NC_Serial.cc:77
void inq_varnatts_impl(const std::string &variable_name, int &result) const
Get the number of attributes of a variable.
Definition: NC_Serial.cc:544
void inq_varid_impl(const std::string &variable_name, bool &exists) const
Finds a variable and sets the "exists" flag.
Definition: NC_Serial.cc:565
void inq_attname_impl(const std::string &variable_name, unsigned int n, std::string &result) const
Gets the name of a numbered attribute.
Definition: NC_Serial.cc:797
virtual void def_var_impl(const std::string &name, io::Type nctype, const std::vector< std::string > &dims) const
Define a variable.
Definition: NC_Serial.cc:252
virtual void set_compression_level_impl(int level) const
Definition: NC_Serial.cc:71
void close_impl()
Close a NetCDF file.
Definition: NC_Serial.cc:109
NC_Serial(MPI_Comm com)
Definition: NC_Serial.cc:55
void del_att_impl(const std::string &variable_name, const std::string &att_name) const
Definition: NC_Serial.cc:894
void inq_nvars_impl(int &result) const
Get the number of variables.
Definition: NC_Serial.cc:473
void sync_impl() const
Definition: NC_Serial.cc:125
void get_var_double(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, const std::vector< unsigned int > &imap, double *ip, bool transposed) const
Get variable data.
Definition: NC_Serial.cc:296
void inq_dimlen_impl(const std::string &dimension_name, unsigned int &result) const
Get a dimension length.
Definition: NC_Serial.cc:201
void redef_impl() const
Enter define mode.
Definition: NC_Serial.cc:156
void inq_dimid_impl(const std::string &dimension_name, bool &exists) const
Definition: NC_Serial.cc:185
virtual ~NC_Serial()
Definition: NC_Serial.cc:60
int get_varid(const std::string &variable_name) const
Definition: NC_Serial.cc:917
void inq_unlimdim_impl(std::string &result) const
Get an unlimited dimension.
Definition: NC_Serial.cc:224
void inq_atttype_impl(const std::string &variable_name, const std::string &att_name, io::Type &result) const
Gets the type of an attribute.
Definition: NC_Serial.cc:825
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
static void check(const ErrorLocation &where, int return_code)
Prints an error message; for debugging.
Definition: NC4_Par.cc:36
static void check_and_abort(MPI_Comm com, const ErrorLocation &where, int return_code)
call MPI_Abort() if a NetCDF call failed
Definition: NC_Serial.cc:48
static void get_att_string(int ncid, int varid, const std::string &att_name, std::string &result)
Definition: NC4File.cc:307
static void get_att_text(int ncid, int varid, const std::string &att_name, std::string &result)
Definition: NC4File.cc:290
static const double k
Definition: exactTestP.cc:42
std::string join(const std::vector< std::string > &strings, const std::string &separator)
Concatenate strings, inserting separator between elements.
static pism::io::Type nc_type_to_pism_type(int input)
static nc_type pism_type_to_nc_type(pism::io::Type input)
int count
Definition: test_cube.c:16