PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
pism_utilities.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2021, 2022 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 
20 #include "pism_utilities.hh"
21 
22 #include <cstdarg> // va_list, va_start(), va_end()
23 #include <sstream> // istringstream, ostringstream
24 #include <cstdio> // vsnprintf
25 #include <cassert> // assert
26 
27 #include <mpi.h> // MPI_Get_library_version
28 #include <fftw3.h> // fftw_version
29 #include <gsl/gsl_version.h> // GSL_VERSION
30 
31 // All recent NetCDF versions depend on HDF5, so we should be able to include hdf5.h to
32 // record the version used by PISM.
33 #include <hdf5.h> // H5_VERS_INFO
34 
35 #include "pism/pism_config.hh" // Pism_USE_XXX, version info
36 
37 // The following is a stupid kludge necessary to make NetCDF 4.x work in
38 // serial mode in an MPI program:
39 #ifndef MPI_INCLUDED
40 #define MPI_INCLUDED 1
41 #endif
42 #include <netcdf.h> // nc_inq_libvers
43 
44 #if (Pism_USE_PROJ==1)
45 #include "pism/util/Proj.hh" // pj_release
46 #endif
47 
48 #if (Pism_USE_JANSSON==1)
49 #include <jansson.h> // JANSSON_VERSION
50 #endif
51 
52 #include <petsctime.h> // PetscTime
53 
54 #include <cstdlib> // strtol(), strtod()
55 
56 #include "error_handling.hh"
57 
58 namespace pism {
59 
60 std::string string_strip(const std::string &input) {
61  if (input.empty()) {
62  return "";
63  }
64 
65  std::string tmp = input;
66 
67  // strip leading spaces
68  tmp.erase(0, tmp.find_first_not_of(" \t"));
69 
70  // strip trailing spaces
71  tmp.substr(tmp.find_last_not_of(" \t"));
72 
73  return tmp;
74 }
75 
76 //! Returns true if `str` ends with `suffix` and false otherwise.
77 bool ends_with(const std::string &str, const std::string &suffix) {
78  if (suffix.size() > str.size()) {
79  return false;
80  }
81 
82  return (str.rfind(suffix) + suffix.size() == str.size());
83 }
84 
85 template <class T>
86 std::string join_impl(const T& input, const std::string& separator) {
87  auto j = input.begin();
88  std::string result = *j;
89  ++j;
90  while (j != input.end()) {
91  result += separator + *j;
92  ++j;
93  }
94  return result;
95 }
96 
97 //! Concatenate `strings`, inserting `separator` between elements.
98 std::string join(const std::vector<std::string> &strings, const std::string &separator) {
99  return join_impl(strings, separator);
100 }
101 
102 std::string set_join(const std::set<std::string> &input, const std::string& separator) {
103  return join_impl(input, separator);
104 }
105 
106 /*!
107  * Replace all occurrences of the character `from` with `to` and return the resulting string.
108  *
109  * We could use std::regex_replace(), but this will do for now.
110  */
111 std::string replace_character(const std::string &input, char from, char to) {
112  std::string result;
113  for (const char &c : input) {
114  result += c == from ? to : c;
115  }
116  return result;
117 }
118 
119 
120 //! Transform a `separator`-separated list (a string) into a vector of strings.
121 std::vector<std::string> split(const std::string &input, char separator) {
122  std::istringstream input_list(input);
123  std::string token;
124  std::vector<std::string> result;
125 
126  while (getline(input_list, token, separator)) {
127  if (not token.empty()) {
128  result.emplace_back(token);
129  }
130  }
131  return result;
132 }
133 
134 //! Transform a `separator`-separated list (a string) into a set of strings.
135 std::set<std::string> set_split(const std::string &input, char separator) {
136  std::set<std::string> result;
137  for (const auto &token : split(input, separator)) {
138  result.insert(token);
139  }
140  return result;
141 }
142 
143 //! Checks if a vector of doubles is strictly increasing.
144 bool is_increasing(const std::vector<double> &a) {
145  int len = (int)a.size();
146  for (int k = 0; k < len-1; k++) {
147  if (a[k] >= a[k+1]) {
148  return false;
149  }
150  }
151  return true;
152 }
153 
154 bool member(const std::string &string, const std::set<std::string> &set) {
155  return (set.find(string) != set.end());
156 }
157 
158 void GlobalReduce(MPI_Comm comm, double *local, double *result, int count, MPI_Op op) {
159  assert(local != result);
160  int err = MPI_Allreduce(local, result, count, MPI_DOUBLE, op, comm);
161  PISM_C_CHK(err, 0, "MPI_Allreduce");
162 }
163 
164 void GlobalReduce(MPI_Comm comm, int *local, int *result, int count, MPI_Op op) {
165  assert(local != result);
166  int err = MPI_Allreduce(local, result, count, MPI_INT, op, comm);
167  PISM_C_CHK(err, 0, "MPI_Allreduce");
168 }
169 
170 void GlobalMin(MPI_Comm comm, double *local, double *result, int count) {
171  GlobalReduce(comm, local, result, count, MPI_MIN);
172 }
173 
174 void GlobalMax(MPI_Comm comm, double *local, double *result, int count) {
175  GlobalReduce(comm, local, result, count, MPI_MAX);
176 }
177 
178 void GlobalMax(MPI_Comm comm, int *local, int *result, int count) {
179  GlobalReduce(comm, local, result, count, MPI_MAX);
180 }
181 
182 void GlobalSum(MPI_Comm comm, double *local, double *result, int count) {
183  GlobalReduce(comm, local, result, count, MPI_SUM);
184 }
185 
186 void GlobalSum(MPI_Comm comm, int *local, int *result, int count) {
187  GlobalReduce(comm, local, result, count, MPI_SUM);
188 }
189 
190 unsigned int GlobalSum(MPI_Comm comm, unsigned int input) {
191  unsigned int result;
192  int err = MPI_Allreduce(&input, &result, 1, MPI_UNSIGNED, MPI_SUM, comm);
193  PISM_C_CHK(err, 0, "MPI_Allreduce");
194  return result;
195 }
196 
197 int GlobalSum(MPI_Comm comm, int input) {
198  int result;
199  int err = MPI_Allreduce(&input, &result, 1, MPI_INT, MPI_SUM, comm);
200  PISM_C_CHK(err, 0, "MPI_Allreduce");
201  return result;
202 }
203 
204 double GlobalMin(MPI_Comm comm, double local) {
205  double result;
206  GlobalMin(comm, &local, &result, 1);
207  return result;
208 }
209 
210 double GlobalMax(MPI_Comm comm, double local) {
211  double result;
212  GlobalMax(comm, &local, &result, 1);
213  return result;
214 }
215 
216 double GlobalSum(MPI_Comm comm, double local) {
217  double result;
218  GlobalSum(comm, &local, &result, 1);
219  return result;
220 }
221 
222 static const int TEMPORARY_STRING_LENGTH = 32768;
223 
224 std::string version() {
225  char buffer[TEMPORARY_STRING_LENGTH];
226  std::string result;
227 
228  result += pism::printf("PISM (%s)\n", pism::revision);
229  result += pism::printf("CMake %s.\n", pism::cmake_version);
230 
231  PetscGetVersion(buffer, TEMPORARY_STRING_LENGTH);
232  result += buffer;
233  result += "\n";
234  result += pism::printf("PETSc configure: %s\n", pism::petsc_configure_flags);
235 
236  // OpenMPI added MPI_Get_library_version in version 1.7 (relatively recently).
237 #ifdef OPEN_MPI
238  result += pism::printf("OpenMPI %d.%d.%d\n",
239  OMPI_MAJOR_VERSION,
240  OMPI_MINOR_VERSION,
241  OMPI_RELEASE_VERSION);
242 #else
243  // Assume that other MPI libraries implement this part of the MPI-3 standard...
244  int string_length = TEMPORARY_STRING_LENGTH;
245  MPI_Get_library_version(buffer, &string_length);
246  result += buffer;
247 #endif
248 
249  result += pism::printf("%s\n", H5_VERS_INFO);
250  result += pism::printf("NetCDF %s.\n", nc_inq_libvers());
251  result += pism::printf("FFTW %s.\n", fftw_version);
252  result += pism::printf("GSL %s.\n", GSL_VERSION);
253 
254 #if (Pism_USE_PROJ==1)
255  result += pism::printf("PROJ %s.\n", pj_release);
256 #endif
257 
258 #if (Pism_USE_JANSSON==1)
259  result += pism::printf("Jansson %s.\n", JANSSON_VERSION);
260 #endif
261 
262 #if (Pism_BUILD_PYTHON_BINDINGS==1)
263  result += pism::printf("SWIG %s.\n", pism::swig_version);
264  result += pism::printf("petsc4py %s.\n", pism::petsc4py_version);
265 #endif
266 
267  return result;
268 }
269 
270 
271 //! Return time since the beginning of the run, in hours.
272 double wall_clock_hours(MPI_Comm com, double start_time) {
273  int rank = 0;
274  double result = 0.0;
275 
276  MPI_Comm_rank(com, &rank);
277 
278  ParallelSection rank0(com);
279  try {
280  if (rank == 0) {
281  const double seconds_per_hour = 3600.0;
282  result = (get_time() - start_time) / seconds_per_hour;
283  }
284  } catch (...) {
285  rank0.failed();
286  }
287  rank0.check();
288 
289  MPI_Bcast(&result, 1, MPI_DOUBLE, 0, com);
290 
291  return result;
292 }
293 
294 //! Creates a time-stamp used for the history NetCDF attribute.
295 std::string timestamp(MPI_Comm com) {
296  time_t now;
297  tm tm_now;
298  char date_str[50];
299  now = time(NULL);
300  localtime_r(&now, &tm_now);
301  // Format specifiers for strftime():
302  // %F = ISO date format, %T = Full 24 hour time, %Z = Time Zone name
303  strftime(date_str, sizeof(date_str), "%F %T %Z", &tm_now);
304 
305  MPI_Bcast(date_str, 50, MPI_CHAR, 0, com);
306 
307  return std::string(date_str);
308 }
309 
310 //! Creates a string with the user name, hostname and the time-stamp (for history strings).
311 std::string username_prefix(MPI_Comm com) {
312  PetscErrorCode ierr;
313 
314  char username[50];
315  ierr = PetscGetUserName(username, sizeof(username));
316  PISM_CHK(ierr, "PetscGetUserName");
317  if (ierr != 0) {
318  username[0] = '\0';
319  }
320  char hostname[100];
321  ierr = PetscGetHostName(hostname, sizeof(hostname));
322  PISM_CHK(ierr, "PetscGetHostName");
323  if (ierr != 0) {
324  hostname[0] = '\0';
325  }
326 
327  std::ostringstream message;
328  message << username << "@" << hostname << " " << timestamp(com) << ": ";
329 
330  std::string result = message.str();
331  unsigned int length = result.size();
332  MPI_Bcast(&length, 1, MPI_UNSIGNED, 0, com);
333 
334  result.resize(length);
335  MPI_Bcast(&result[0], length, MPI_CHAR, 0, com);
336 
337  return result;
338 }
339 
340 //! \brief Uses argc and argv to create the string with current PISM
341 //! command-line arguments.
342 std::string args_string() {
343  int argc;
344  char **argv;
345  PetscErrorCode ierr = PetscGetArgs(&argc, &argv);
346  PISM_CHK(ierr, "PetscGetArgs");
347 
348  std::string cmdstr;
349  std::string argument;
350  for (int j = 0; j < argc; j++) {
351  argument = argv[j];
352 
353  // enclose arguments containing spaces with double quotes:
354  if (argument.find(' ') != std::string::npos) {
355  argument = "\"" + argument + "\"";
356  }
357 
358  cmdstr += std::string(" ") + argument;
359  }
360  cmdstr += "\n";
361 
362  return cmdstr;
363 }
364 
365 //! \brief Adds a suffix to a filename.
366 /*!
367  * Returns filename + separator + suffix + .nc if the original filename had the
368  * .nc suffix, otherwise filename + separator. If the old filename had the form
369  * "name + separator + more stuff + .nc", then removes the string after the
370  * separator.
371  */
372 std::string filename_add_suffix(const std::string &filename,
373  const std::string &separator,
374  const std::string &suffix) {
375  std::string basename = filename;
376  std::string result;
377 
378  // find where the separator begins:
379  std::string::size_type j = basename.rfind(separator);
380  if (j == std::string::npos) {
381  j = basename.rfind(".nc");
382  }
383 
384  // if the separator was not found, find the .nc suffix:
385  if (j == std::string::npos) {
386  j = basename.size();
387  }
388 
389  // cut off everything starting from the separator (or the .nc suffix):
390  basename.resize(static_cast<int>(j));
391 
392  result = basename + separator + suffix;
393 
394  if (ends_with(filename, ".nc")) {
395  result += ".nc";
396  }
397 
398  return result;
399 }
400 
401 double get_time() {
402  PetscLogDouble result;
403  PetscErrorCode ierr = PetscTime(&result); PISM_CHK(ierr, "PetscTime");
404  return result;
405 }
406 
407 std::string printf(const char *format, ...) {
408  std::string result(1024, ' ');
409  va_list arglist;
410  size_t length;
411 
412  va_start(arglist, format);
413  if((length = vsnprintf(&result[0], result.size(), format, arglist)) > result.size()) {
414  result.reserve(length);
415  vsnprintf(&result[0], result.size(), format, arglist);
416  }
417  va_end(arglist);
418  return result.substr(0, length);
419 }
420 
421 /*!
422  * Validate a format string. In this application a format string should contain `%` exactly
423  * once, followed by `s` (i.e. `%s`).
424  *
425  * Throws RuntimeError if the provided string is invalid.
426  */
427 void validate_format_string(const std::string &format) {
428  if (format.find("%s") == std::string::npos) {
429  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "format string %s does not contain %%s",
430  format.c_str());
431  }
432 
433  if (format.find('%') != format.rfind('%')) {
434  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "format string %s contains more than one %%",
435  format.c_str());
436  }
437 }
438 
439 double vector_min(const std::vector<double> &input) {
440  double my_min = input[0];
441  for (auto x : input) {
442  my_min = std::min(x, my_min);
443  }
444  return my_min;
445 }
446 
447 double vector_max(const std::vector<double> &input) {
448  double my_max = input[0];
449  for (auto x : input) {
450  my_max = std::max(x, my_max);
451  }
452  return my_max;
453 }
454 
455 
456 /*!
457  * Fletcher's checksum
458  *
459  * See https://en.wikipedia.org/wiki/Fletcher%27s_checksum#Optimizations
460  */
461 uint64_t fletcher64(const uint32_t *data, size_t length) {
462  // Accumulating a sum of block_size unsigned 32-bit integers in an unsigned 64-bit
463  // integer will not lead to an overflow.
464  //
465  // This constant is found by solving n * (n + 1) / 2 * (2^32 - 1) < (2^64 - 1).
466  static const size_t block_size = 92681;
467 
468  uint64_t c0 = 0;
469  uint64_t c1 = 0;
470  while (length != 0) {
471  size_t block = std::min(block_size, length);
472 
473  for (size_t i = 0; i < block; ++i) {
474  c0 = c0 + *data++;
475  c1 = c1 + c0;
476  }
477 
478  c0 = c0 % UINT32_MAX;
479  c1 = c1 % UINT32_MAX;
480 
481  length = length > block_size ? length - block_size : 0;
482  }
483  return (c1 << 32 | c0);
484 }
485 
486 /*!
487  * Call fletcher64() to compute the checksum and print it.
488  */
489 void print_checksum(MPI_Comm com,
490  const std::vector<double> &data,
491  const char *label) {
492  int rank = 0;
493  MPI_Comm_rank(com, &rank);
494 
495  uint64_t sum = fletcher64((uint32_t*)data.data(), data.size() * 2);
496 
497  PetscPrintf(PETSC_COMM_SELF, "[%d] %s: %016llx\n",
498  rank, label, (unsigned long long int)sum);
499 }
500 
501 void print_vector(MPI_Comm com,
502  const std::vector<double> &data,
503  const char *label) {
504  int rank = 0;
505  MPI_Comm_rank(com, &rank);
506 
507  std::vector<std::string> tmp;
508  for (const auto &f : data) {
509  tmp.emplace_back(pism::printf("%f", f));
510  }
511 
512  auto str = join(tmp, ",");
513 
514  PetscPrintf(PETSC_COMM_SELF, "[%d] %s: %s\n",
515  rank, label, str.c_str());
516 }
517 
518 void print_vector(MPI_Comm com,
519  const std::vector<int> &data,
520  const char *label) {
521  int rank = 0;
522  MPI_Comm_rank(com, &rank);
523 
524  std::vector<std::string> tmp;
525  for (const auto &f : data) {
526  tmp.emplace_back(pism::printf("%d", f));
527  }
528 
529  auto str = join(tmp, ",");
530 
531  PetscPrintf(PETSC_COMM_SELF, "[%d] %s: %s\n",
532  rank, label, str.c_str());
533 }
534 
535 /*!
536  * Compute water column pressure vertically-averaged over the height of an ice cliff at a
537  * margin.
538  */
539 double average_water_column_pressure(double ice_thickness, double bed,
540  double floatation_level, double rho_ice,
541  double rho_water, double g) {
542 
543  double ice_bottom = std::max(bed, floatation_level - rho_ice / rho_water * ice_thickness);
544  double water_column_height = std::max(floatation_level - ice_bottom, 0.0);
545 
546  if (ice_thickness > 0.0) {
547  return 0.5 * rho_water * g * pow(water_column_height, 2.0) / ice_thickness;
548  }
549  return 0.0;
550 }
551 
552 double parse_number(const std::string &input) {
553  char *endptr = NULL;
554  double result = strtod(input.c_str(), &endptr);
555  if (*endptr != '\0') {
557  "Can't parse %s (expected a floating point number)",
558  input.c_str());
559  }
560  return result;
561 }
562 
563 long int parse_integer(const std::string &input) {
564  char *endptr = NULL;
565  long int result = strtol(input.c_str(), &endptr, 10);
566  if (*endptr != '\0') {
568  "Can't parse %s (expected an integer)",
569  input.c_str());
570  }
571  return result;
572 }
573 
574 } // end of namespace pism
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
#define PISM_C_CHK(errcode, success, name)
const double rho_ice
Definition: exactTestK.c:31
double get_time()
bool is_increasing(const std::vector< double > &a)
Checks if a vector of doubles is strictly increasing.
double parse_number(const std::string &input)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
double average_water_column_pressure(double ice_thickness, double bed, double floatation_level, double rho_ice, double rho_water, double g)
double wall_clock_hours(MPI_Comm com, double start_time)
Return time since the beginning of the run, in hours.
static double start_time(const Config &config, const Logger &log, const File *file, const std::string &reference_date, const std::string &calendar, const units::Unit &time_units)
Definition: Time.cc:347
void print_checksum(MPI_Comm com, const std::vector< double > &data, const char *label)
static const double g
Definition: exactTestP.cc:39
bool ends_with(const std::string &str, const std::string &suffix)
Returns true if str ends with suffix and false otherwise.
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
std::string set_join(const std::set< std::string > &input, const std::string &separator)
std::string printf(const char *format,...)
void print_vector(MPI_Comm com, const std::vector< double > &data, const char *label)
uint64_t fletcher64(const uint32_t *data, size_t length)
static const double k
Definition: exactTestP.cc:45
void validate_format_string(const std::string &format)
std::set< std::string > set_split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a set of strings.
double vector_min(const std::vector< double > &input)
std::string string_strip(const std::string &input)
std::string filename_add_suffix(const std::string &filename, const std::string &separator, const std::string &suffix)
Adds a suffix to a filename.
std::string version()
std::string timestamp(MPI_Comm com)
Creates a time-stamp used for the history NetCDF attribute.
bool member(const std::string &string, const std::set< std::string > &set)
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double c1
Definition: exactTestP.cc:47
void GlobalMin(MPI_Comm comm, double *local, double *result, int count)
std::string args_string()
Uses argc and argv to create the string with current PISM command-line arguments.
long int parse_integer(const std::string &input)
static const int TEMPORARY_STRING_LENGTH
std::string join(const std::vector< std::string > &strings, const std::string &separator)
Concatenate strings, inserting separator between elements.
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
std::string username_prefix(MPI_Comm com)
Creates a string with the user name, hostname and the time-stamp (for history strings).
std::string join_impl(const T &input, const std::string &separator)
void GlobalReduce(MPI_Comm comm, double *local, double *result, int count, MPI_Op op)
std::string replace_character(const std::string &input, char from, char to)
double sum(const IceModelVec2S &input)
Sums up all the values in an IceModelVec2S object. Ignores ghosts.
double vector_max(const std::vector< double > &input)
std::vector< std::string > split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a vector of strings.
int count
Definition: test_cube.c:16