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