PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
iceModelVec.cc
Go to the documentation of this file.
1 // Copyright (C) 2008--2021 Ed Bueler, Constantine Khroulev, and David Maxwell
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 <cassert>
20 
21 #include <petscdraw.h>
22 
23 #include "iceModelVec.hh"
24 #include "pism/util/IceModelVec2V.hh"
25 #include "pism/util/IceModelVec_impl.hh"
26 
27 #include "Time.hh"
28 #include "IceGrid.hh"
29 #include "ConfigInterface.hh"
30 
31 #include "error_handling.hh"
32 #include "io/io_helpers.hh"
33 #include "pism/util/Logger.hh"
34 #include "pism/util/Profiling.hh"
35 #include "pism/util/petscwrappers/VecScatter.hh"
36 #include "pism/util/petscwrappers/Viewer.hh"
37 #include "pism/util/Mask.hh"
38 #include "pism/util/IceModelVec2CellType.hh"
39 #include "pism/util/Context.hh"
40 #include "pism/util/VariableMetadata.hh"
41 #include "pism/util/io/File.hh"
42 #include "pism_utilities.hh"
43 
44 namespace pism {
45 
46 static void global_to_local(petsc::DM &dm, Vec source, Vec destination) {
47  PetscErrorCode ierr;
48 
49  ierr = DMGlobalToLocalBegin(dm, source, INSERT_VALUES, destination);
50  PISM_CHK(ierr, "DMGlobalToLocalBegin");
51 
52  ierr = DMGlobalToLocalEnd(dm, source, INSERT_VALUES, destination);
53  PISM_CHK(ierr, "DMGlobalToLocalEnd");
54 }
55 
57  const std::string &name,
58  IceModelVecKind ghostedp,
59  size_t dof,
60  size_t stencil_width,
61  const std::vector<double> &zlevels) {
62  m_impl = new Impl();
63  m_array = nullptr;
64 
65  m_impl->name = name;
66  m_impl->grid = grid;
67  m_impl->ghosted = (ghostedp == WITH_GHOSTS);
68  m_impl->dof = dof;
69  m_impl->zlevels = zlevels;
70 
71  auto max_stencil_width = grid->ctx()->config()->get_number("grid.max_stencil_width");
72  if ((dof != 1) or (stencil_width > max_stencil_width)) {
73  // use the requested stencil width *if* we have to
75  } else {
76  // otherwise use the "standard" stencil width
77  m_impl->da_stencil_width = max_stencil_width;
78  }
79 
80  auto system = m_impl->grid->ctx()->unit_system();
81  if (m_impl->dof > 1) {
82  // dof > 1: this is a 2D vector
83  using pism::printf;
84  for (unsigned int j = 0; j < m_impl->dof; ++j) {
85  m_impl->metadata.push_back({system, printf("%s[%d]", name.c_str(), j)});
86  }
87  } else {
88  // both 2D and 3D vectors
89  m_impl->metadata = {{system, name, zlevels}};
90  }
91 
92  if (zlevels.size() > 1) {
93  m_impl->bsearch_accel = gsl_interp_accel_alloc();
94  if (m_impl->bsearch_accel == NULL) {
96  "Failed to allocate a GSL interpolation accelerator");
97  }
98  }
99 }
100 
102  assert(m_impl->access_counter == 0);
103 
104  if (m_impl->bsearch_accel != nullptr) {
105  gsl_interp_accel_free(m_impl->bsearch_accel);
106  m_impl->bsearch_accel = nullptr;
107  }
108 
109  delete m_impl;
110  m_impl = nullptr;
111 }
112 
113 //! \brief Get the object state counter.
114 /*!
115  * This method returns the "revision number" of an IceModelVec.
116  *
117  * It can be used to determine it a field was updated and if a certain
118  * computation needs to be re-done. One example is computing the smoothed bed
119  * for the SIA computation, which is only necessary if the bed deformation code
120  * fired.
121  *
122  * See also inc_state_counter().
123  */
125  return m_impl->state_counter;
126 }
127 
129  return m_impl->grid;
130 }
131 
132 unsigned int IceModelVec::ndof() const {
133  return m_impl->dof;
134 }
135 
136 std::vector<double> IceModelVec::levels() const {
137  return m_impl->zlevels;
138 }
139 
140 //! \brief Increment the object state counter.
141 /*!
142  * See the documentation of get_state_counter(). This method is the
143  * *only* way to manually increment the state counter. It is also
144  * automatically updated by IceModelVec methods that are known to
145  * change stored values.
146  */
149 }
150 
151 //! Returns the number of spatial dimensions.
152 unsigned int IceModelVec::ndims() const {
153  return m_impl->zlevels.size() > 1 ? 3 : 2;
154 }
155 
156 std::vector<int> IceModelVec::shape() const {
157 
158  auto grid = m_impl->grid;
159 
160  if (ndims() == 3) {
161  return {(int)grid->My(), (int)grid->Mx(), (int)levels().size()};
162  }
163 
164  if (ndof() == 1) {
165  return {(int)grid->My(), (int)grid->Mx()};
166  }
167 
168  return {(int)grid->My(), (int)grid->Mx(), (int)ndof()};
169 }
170 
171 //! Set the time independent flag for all variables corresponding to this IceModelVec instance.
172 /** A "time independent" IceModelVec will be saved to a NetCDF
173  variable which does not depend on the "time" dimension.
174  */
176  for (unsigned int j = 0; j < m_impl->dof; ++j) {
177  m_impl->metadata[j].set_time_independent(flag);
178  }
179 }
180 
183 }
184 
185 
186 //! Result: min <- min(v[j]), max <- max(v[j]).
187 /*!
188 PETSc manual correctly says "VecMin and VecMax are collective on Vec" but
189 GlobalMax,GlobalMin \e are needed, when m_impl->ghosted==true, to get correct
190 values because Vecs created with DACreateLocalVector() are of type
191 VECSEQ and not VECMPI. See src/trypetsc/localVecMax.c.
192  */
193 std::array<double,2> IceModelVec::range() const {
194  PetscErrorCode ierr;
195 
196  double min{0.0};
197  ierr = VecMin(vec(), NULL, &min);
198  PISM_CHK(ierr, "VecMin");
199 
200  double max{0.0};
201  ierr = VecMax(vec(), NULL, &max);
202  PISM_CHK(ierr, "VecMax");
203 
204  if (m_impl->ghosted) {
205  // needs a reduce operation; use GlobalMin and GlobalMax;
206  min = GlobalMin(m_impl->grid->com, min);
207  max = GlobalMax(m_impl->grid->com, max);
208  }
209  return {min, max};
210 }
211 
212 /** Convert from `int` to PETSc's `NormType`.
213  *
214  *
215  * @param[in] input norm type as an integer
216  *
217  * @return norm type as PETSc's `NormType`.
218  */
219 static NormType int_to_normtype(int input) {
220  switch (input) {
221  case NORM_1:
222  return NORM_1;
223  case NORM_2:
224  return NORM_2;
225  case NORM_INFINITY:
226  return NORM_INFINITY;
227  default:
228  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid norm type: %d",
229  input);
230  }
231 }
232 
233 //! Result: v <- v + alpha * x. Calls VecAXPY.
234 void IceModelVec::add(double alpha, const IceModelVec &x) {
235  checkCompatibility("add", x);
236 
237  PetscErrorCode ierr = VecAXPY(vec(), alpha, x.vec());
238  PISM_CHK(ierr, "VecAXPY");
239 
240  inc_state_counter(); // mark as modified
241 }
242 
243 //! Result: v[j] <- v[j] + alpha for all j. Calls VecShift.
244 void IceModelVec::shift(double alpha) {
245  PetscErrorCode ierr = VecShift(vec(), alpha);
246  PISM_CHK(ierr, "VecShift");
247 
248  inc_state_counter(); // mark as modified
249 }
250 
251 //! Result: v <- v * alpha. Calls VecScale.
252 void IceModelVec::scale(double alpha) {
253  PetscErrorCode ierr = VecScale(vec(), alpha);
254  PISM_CHK(ierr, "VecScale");
255 
256  inc_state_counter(); // mark as modified
257 }
258 
259 //! Copies v to a global vector 'destination'. Ghost points are discarded.
260 /*! This is potentially dangerous: make sure that `destination` has the same
261  dimensions as the current IceModelVec.
262 
263  DMLocalToGlobalBegin/End is broken in PETSc 3.5, so we roll our
264  own.
265  */
266 void IceModelVec::copy_to_vec(std::shared_ptr<petsc::DM> destination_da,
267  petsc::Vec &destination) const {
268  // m_dof > 1 for vector, staggered grid 2D fields, etc. In this case
269  // zlevels.size() == 1. For 3D fields, m_dof == 1 (all 3D fields are
270  // scalar) and zlevels.size() corresponds to dof of the underlying PETSc
271  // DM object. So we want the bigger of the two numbers here.
272  unsigned int N = std::max((size_t)m_impl->dof, m_impl->zlevels.size());
273 
274  this->get_dof(destination_da, destination, 0, N);
275 }
276 
277 void IceModelVec::get_dof(std::shared_ptr<petsc::DM> da_result,
278  petsc::Vec &result,
279  unsigned int start, unsigned int count) const {
280  if (start >= m_impl->dof) {
281  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid argument (start); got %d", start);
282  }
283 
284  petsc::DMDAVecArrayDOF tmp_res(da_result, result), tmp_v(dm(), vec());
285 
286  double
287  ***result_a = static_cast<double***>(tmp_res.get()),
288  ***source_a = static_cast<double***>(tmp_v.get());
289 
290  ParallelSection loop(m_impl->grid->com);
291  try {
292  for (Points p(*m_impl->grid); p; p.next()) {
293  const int i = p.i(), j = p.j();
294  PetscErrorCode ierr = PetscMemcpy(result_a[j][i], &source_a[j][i][start],
295  count*sizeof(PetscScalar));
296  PISM_CHK(ierr, "PetscMemcpy");
297  }
298  } catch (...) {
299  loop.failed();
300  }
301  loop.check();
302 }
303 
304 void IceModelVec::set_dof(std::shared_ptr<petsc::DM> da_source, petsc::Vec &source,
305  unsigned int start, unsigned int count) {
306  if (start >= m_impl->dof) {
307  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid argument (start); got %d", start);
308  }
309 
310  petsc::DMDAVecArrayDOF tmp_src(da_source, source), tmp_v(dm(), vec());
311 
312  double
313  ***source_a = static_cast<double***>(tmp_src.get()),
314  ***result_a = static_cast<double***>(tmp_v.get());
315 
316  ParallelSection loop(m_impl->grid->com);
317  try {
318  for (Points p(*m_impl->grid); p; p.next()) {
319  const int i = p.i(), j = p.j();
320  PetscErrorCode ierr = PetscMemcpy(&result_a[j][i][start], source_a[j][i],
321  count*sizeof(PetscScalar));
322  PISM_CHK(ierr, "PetscMemcpy");
323  }
324  } catch (...) {
325  loop.failed();
326  }
327  loop.check();
328 
329  inc_state_counter(); // mark as modified
330 }
331 
332 //! @brief Get the stencil width of the current IceModelVec. Returns 0
333 //! if ghosts are not available.
334 unsigned int IceModelVec::stencil_width() const {
335  if (m_impl->ghosted) {
336  return m_impl->da_stencil_width;
337  }
338 
339  return 0;
340 }
341 
343  if (m_impl->v.get() == nullptr) {
344  PetscErrorCode ierr = 0;
345  if (m_impl->ghosted) {
346  ierr = DMCreateLocalVector(*dm(), m_impl->v.rawptr());
347  PISM_CHK(ierr, "DMCreateLocalVector");
348  } else {
349  ierr = DMCreateGlobalVector(*dm(), m_impl->v.rawptr());
350  PISM_CHK(ierr, "DMCreateGlobalVector");
351  }
352  }
353  return m_impl->v;
354 }
355 
356 std::shared_ptr<petsc::DM> IceModelVec::dm() const {
357  if (m_impl->da == nullptr) {
358  // dof > 1 for vector, staggered grid 2D fields, etc. In this case zlevels.size() ==
359  // 1. For 3D fields, dof == 1 (all 3D fields are scalar) and zlevels.size()
360  // corresponds to dof of the underlying PETSc DM object.
361  auto da_dof = std::max(m_impl->zlevels.size(), (size_t)m_impl->dof);
362 
363  // initialize the da member:
364  m_impl->da = grid()->get_dm(da_dof, m_impl->da_stencil_width);
365  }
366  return m_impl->da;
367 }
368 
369 //! Sets the variable name to `name`.
370 /**
371  * This is the "overall" name of a field. This is **not** the same as the
372  * NetCDF variable name. Use `metadata(...).set_name(...)` to set that.
373  */
374 void IceModelVec::set_name(const std::string &name) {
375  m_impl->name = name;
376 }
377 
378 //! @brief Get the name of an IceModelVec object.
379 /**
380  * This is the name used to refer to this object in PISM (e.g. via the
381  * Vars class), **not** the name of the corresponding NetCDF variable.
382  * (The problem is that one IceModelVec instance may correspond to
383  * several NetCDF variables. Use metadata(...).get_name() to get the
384  * name of NetCDF variables an IceModelVec is saved to.)
385  */
386 const std::string& IceModelVec::get_name() const {
387  return m_impl->name;
388 }
389 
390 //! Sets NetCDF attributes of an IceModelVec object.
391 /*! Call set_attrs("new pism_intent", "new long name", "new units", "") if a
392  variable does not have a standard name. Similarly, by putting "" in an
393  appropriate spot, it is possible tp leave long_name, units or pism_intent
394  unmodified.
395 
396  If units != "", this also resets glaciological_units, so that they match
397  internal units.
398  */
399 void IceModelVec::set_attrs(const std::string &pism_intent,
400  const std::string &long_name,
401  const std::string &units,
402  const std::string &glaciological_units,
403  const std::string &standard_name,
404  unsigned int component) {
405 
406  metadata(component)["long_name"] = long_name;
407 
408  metadata(component)["units"] = units;
409 
410  if (not m_impl->grid->ctx()->config()->get_flag("output.use_MKS")) {
411  metadata(component)["glaciological_units"] = glaciological_units;
412  }
413 
414  metadata(component)["pism_intent"] = pism_intent;
415 
416  metadata(component)["standard_name"] = standard_name;
417 }
418 
419 //! Gets an IceModelVec from a file `file`, interpolating onto the current grid.
420 /*! Stops if the variable was not found and `critical` == true.
421  */
423  double default_value) {
424 
425  bool allow_extrapolation = grid()->ctx()->config()->get_flag("grid.allow_extrapolation");
426 
427  if (ndof() == 1) {
428  if (m_impl->ghosted) {
430  petsc::VecArray tmp_array(tmp);
431 
432  io::regrid_spatial_variable(metadata(0), *grid(), file, flag,
433  m_impl->report_range, allow_extrapolation,
434  default_value, m_impl->interpolation_type,
435  tmp_array.get());
436 
437  global_to_local(*dm(), tmp, vec());
438  } else {
439  petsc::VecArray v_array(vec());
440  io::regrid_spatial_variable(metadata(0), *grid(), file, flag,
441  m_impl->report_range, allow_extrapolation,
442  default_value, m_impl->interpolation_type,
443  v_array.get());
444  }
445  return;
446  }
447 
448  // Get the dof=1, stencil_width=0 DMDA (components are always scalar
449  // and we just need a global Vec):
450  auto da2 = grid()->get_dm(1, 0);
451 
452  // a temporary one-component vector, distributed across processors
453  // the same way v is
454  petsc::TemporaryGlobalVec tmp(da2);
455 
456  for (unsigned int j = 0; j < ndof(); ++j) {
457  {
458  petsc::VecArray tmp_array(tmp);
459  io::regrid_spatial_variable(metadata(j), *grid(), file, flag,
460  m_impl->report_range, allow_extrapolation,
461  default_value, m_impl->interpolation_type,
462  tmp_array.get());
463  }
464 
465  set_dof(da2, tmp, j);
466  }
467 
468  // The calls above only set the values owned by a processor, so we need to
469  // communicate if m_has_ghosts == true:
470  if (m_impl->ghosted) {
471  update_ghosts();
472  }
473 }
474 
475 //! Reads appropriate NetCDF variable(s) into an IceModelVec.
476 void IceModelVec::read_impl(const File &file, const unsigned int time) {
477 
478  Logger::ConstPtr log = grid()->ctx()->log();
479  log->message(4, " Reading %s...\n", m_impl->name.c_str());
480 
481  if (ndof() == 1) {
482  // This takes are of scalar variables (both 2D and 3D).
483  if (m_impl->ghosted) {
485 
486  petsc::VecArray tmp_array(tmp);
487  io::read_spatial_variable(metadata(0), *grid(), file, time, tmp_array.get());
488 
489  global_to_local(*dm(), tmp, vec());
490  } else {
491  petsc::VecArray v_array(vec());
492  io::read_spatial_variable(metadata(0), *grid(), file, time, v_array.get());
493  }
494  return;
495  }
496 
497  // Get the dof=1, stencil_width=0 DMDA (components are always scalar
498  // and we just need a global Vec):
499  auto da2 = grid()->get_dm(1, 0);
500 
501  // A temporary one-component vector, distributed across processors
502  // the same way v is
503  petsc::TemporaryGlobalVec tmp(da2);
504 
505  for (unsigned int j = 0; j < ndof(); ++j) {
506 
507  {
508  petsc::VecArray tmp_array(tmp);
509  io::read_spatial_variable(metadata(j), *grid(), file, time, tmp_array.get());
510  }
511 
512  set_dof(da2, tmp, j);
513  }
514 
515  // The calls above only set the values owned by a processor, so we need to
516  // communicate if m_impl->ghosted is true:
517  if (m_impl->ghosted) {
518  update_ghosts();
519  }
520 }
521 
522 //! \brief Define variables corresponding to an IceModelVec in a file opened using `file`.
523 void IceModelVec::define(const File &file, IO_Type default_type) const {
524  for (unsigned int j = 0; j < ndof(); ++j) {
525  IO_Type type = metadata(j).get_output_type();
526  type = type == PISM_NAT ? default_type : type;
527  io::define_spatial_variable(metadata(j), *m_impl->grid, file, type);
528  }
529 }
530 
531 //! @brief Returns a reference to the SpatialVariableMetadata object
532 //! containing metadata for the compoment N.
534  assert(N < m_impl->dof);
535  return m_impl->metadata[N];
536 }
537 
538 const SpatialVariableMetadata& IceModelVec::metadata(unsigned int N) const {
539  assert(N < m_impl->dof);
540  return m_impl->metadata[N];
541 }
542 
543 //! Writes an IceModelVec to a NetCDF file.
544 void IceModelVec::write_impl(const File &file) const {
545  Logger::ConstPtr log = m_impl->grid->ctx()->log();
546 
547  log->message(4, " Writing %s...\n", m_impl->name.c_str());
548 
549  // The simplest case:
550  if (ndof() == 1) {
551  if (m_impl->ghosted) {
553 
554  this->copy_to_vec(dm(), tmp);
555 
556  petsc::VecArray tmp_array(tmp);
557 
558  io::write_spatial_variable(metadata(0), *grid(), file, tmp_array.get());
559  } else {
560  petsc::VecArray v_array(vec());
561  io::write_spatial_variable(metadata(0), *grid(), file, v_array.get());
562  }
563  return;
564  }
565 
566  // Get the dof=1, stencil_width=0 DMDA (components are always scalar
567  // and we just need a global Vec):
568  auto da2 = m_impl->grid->get_dm(1, 0);
569 
570  // a temporary one-component vector, distributed across processors
571  // the same way v is
572  petsc::TemporaryGlobalVec tmp(da2);
573 
574  for (unsigned int j = 0; j < ndof(); ++j) {
575  get_dof(da2, tmp, j);
576 
577  petsc::VecArray tmp_array(tmp);
578  io::write_spatial_variable(metadata(j), *grid(), file, tmp_array.get());
579  }
580 }
581 
582 //! Dumps a variable to a file, overwriting this file's contents (for debugging).
583 void IceModelVec::dump(const char filename[]) const {
584  File file(m_impl->grid->com, filename,
585  string_to_backend(m_impl->grid->ctx()->config()->get_string("output.format")),
587  m_impl->grid->ctx()->pio_iosys_id());
588 
589  if (not m_impl->metadata[0].get_time_independent()) {
590  io::define_time(file, *m_impl->grid->ctx());
591  io::append_time(file, *m_impl->grid->ctx()->config(), m_impl->grid->ctx()->time()->current());
592  }
593 
594  define(file, PISM_DOUBLE);
595  write(file);
596 }
597 
598 //! Checks if two IceModelVecs have compatible sizes, dimensions and numbers of degrees of freedom.
599 void IceModelVec::checkCompatibility(const char* func, const IceModelVec &other) const {
600  PetscErrorCode ierr;
601  // We have to use PetscInt because of VecGetSizes below.
602  PetscInt X_size, Y_size;
603 
604  if (m_impl->dof != other.m_impl->dof) {
605  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IceModelVec::%s(...): operands have different numbers of degrees of freedom",
606  func);
607  }
608 
609  ierr = VecGetSize(vec(), &X_size);
610  PISM_CHK(ierr, "VecGetSize");
611 
612  ierr = VecGetSize(other.vec(), &Y_size);
613  PISM_CHK(ierr, "VecGetSize");
614 
615  if (X_size != Y_size) {
616  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IceModelVec::%s(...): incompatible Vec sizes (called as %s.%s(%s))",
617  func, m_impl->name.c_str(), func, other.m_impl->name.c_str());
618  }
619 }
620 
621 //! Checks if an IceModelVec is allocated and calls DAVecGetArray.
623 
624  if (m_impl->access_counter < 0) {
625  throw RuntimeError(PISM_ERROR_LOCATION, "IceModelVec::begin_access(): m_access_counter < 0");
626  }
627 
628  if (m_impl->access_counter == 0) {
629  PetscErrorCode ierr;
631  ierr = DMDAVecGetArrayDOF(*dm(), vec(), &m_array);
632  PISM_CHK(ierr, "DMDAVecGetArrayDOF");
633  } else {
634  ierr = DMDAVecGetArray(*dm(), vec(), &m_array);
635  PISM_CHK(ierr, "DMDAVecGetArray");
636  }
637  }
638 
640 }
641 
642 //! Checks if an IceModelVec is allocated and calls DAVecRestoreArray.
644  PetscErrorCode ierr;
645 
646  if (m_array == NULL) {
648  "IceModelVec::end_access(): a == NULL (looks like begin_acces() was not called)");
649  }
650 
651  if (m_impl->access_counter < 0) {
652  throw RuntimeError(PISM_ERROR_LOCATION, "IceModelVec::end_access(): m_access_counter < 0");
653  }
654 
656  if (m_impl->access_counter == 0) {
658  ierr = DMDAVecRestoreArrayDOF(*dm(), vec(), &m_array);
659  PISM_CHK(ierr, "DMDAVecRestoreArrayDOF");
660  } else {
661  ierr = DMDAVecRestoreArray(*dm(), vec(), &m_array);
662  PISM_CHK(ierr, "DMDAVecRestoreArray");
663  }
664  m_array = NULL;
665  }
666 }
667 
668 //! Updates ghost points.
670  PetscErrorCode ierr;
671  if (not m_impl->ghosted) {
672  return;
673  }
674 
675  ierr = DMLocalToLocalBegin(*dm(), vec(), INSERT_VALUES, vec());
676  PISM_CHK(ierr, "DMLocalToLocalBegin");
677 
678  ierr = DMLocalToLocalEnd(*dm(), vec(), INSERT_VALUES, vec());
679  PISM_CHK(ierr, "DMLocalToLocalEnd");
680 }
681 
682 //! Result: v[j] <- c for all j.
683 void IceModelVec::set(const double c) {
684  PetscErrorCode ierr = VecSet(vec(),c);
685  PISM_CHK(ierr, "VecSet");
686 
687  inc_state_counter(); // mark as modified
688 }
689 
690 void IceModelVec::check_array_indices(int i, int j, unsigned int k) const {
691  double ghost_width = 0;
692  if (m_impl->ghosted) {
693  ghost_width = m_impl->da_stencil_width;
694  }
695  // m_impl->dof > 1 for vector, staggered grid 2D fields, etc. In this case
696  // zlevels.size() == 1. For 3D fields, m_impl->dof == 1 (all 3D fields are
697  // scalar) and zlevels.size() corresponds to dof of the underlying PETSc
698  // DM object. So we want the bigger of the two numbers here.
699  unsigned int N = std::max((size_t)m_impl->dof, m_impl->zlevels.size());
700 
701  bool out_of_range = (i < m_impl->grid->xs() - ghost_width) ||
702  (i > m_impl->grid->xs() + m_impl->grid->xm() + ghost_width) ||
703  (j < m_impl->grid->ys() - ghost_width) ||
704  (j > m_impl->grid->ys() + m_impl->grid->ym() + ghost_width) ||
705  (k >= N);
706 
707  if (out_of_range) {
708  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "%s(%d, %d, %d) is out of bounds",
709  m_impl->name.c_str(), i, j, k);
710  }
711 
712  if (m_array == NULL) {
713  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "%s: begin_access() was not called", m_impl->name.c_str());
714  }
715 }
716 
717 namespace vec {
718 namespace details {
719 //! \brief Compute parameters for 2D loop computations involving 3
720 //! IceModelVecs.
721 /*!
722  * Here we assume that z is updated using a local (point-wise) computation
723  * involving x and y.
724  *
725  * "ghosts" is the width of the stencil that can be updated locally.
726  * "scatter" is false if all ghosts can be updated locally.
727  */
728 void compute_params(const IceModelVec* const x, const IceModelVec* const y,
729  const IceModelVec* const z, int &ghosts, bool &scatter) {
730 
731  // We have 2^3=8 cases here (x,y,z having or not having ghosts).
732  if (z->stencil_width() == 0) {
733  // z has no ghosts; we can update everything locally
734  // (This covers 4 cases.)
735  ghosts = 0;
736  scatter = false;
737  } else if (x->stencil_width() == 0 ||
738  y->stencil_width() == 0) {
739  // z has ghosts, but at least one of x and y does not. we have to scatter
740  // ghosts.
741  // (This covers 3 cases.)
742  ghosts = 0;
743  scatter = true;
744  } else {
745  // all of x, y, z have ghosts
746  // (The remaining 8-th case.)
747  if (z->stencil_width() <= x->stencil_width() &&
748  z->stencil_width() <= y->stencil_width()) {
749  // x and y have enough ghosts to update ghosts of z locally
750  ghosts = z->stencil_width();
751  scatter = false;
752  } else {
753  // z has ghosts, but at least one of x and y doesn't have a wide enough
754  // stencil
755  ghosts = 0;
756  scatter = true;
757  }
758  }
759 }
760 
761 } // end of namespace details
762 } // end of namespace vec
763 
764 //! Computes the norm of all the components of an IceModelVec.
765 /*!
766 See comment for range(); because local Vecs are VECSEQ, needs a reduce operation.
767 See src/trypetsc/localVecMax.c.
768  */
769 std::vector<double> IceModelVec::norm(int n) const {
770  std::vector<double> result(m_impl->dof);
771 
772  NormType type = int_to_normtype(n);
773 
774  if (m_impl->dof > 1) {
775  PetscErrorCode ierr = VecStrideNormAll(vec(), type, &result[0]);
776  PISM_CHK(ierr, "VecStrideNormAll");
777  } else {
778  PetscErrorCode ierr = VecNorm(vec(), type, &result[0]);
779  PISM_CHK(ierr, "VecNorm");
780  }
781 
782  if (m_impl->ghosted) {
783  // needs a reduce operation; use GlobalMax() if NORM_INFINITY,
784  // otherwise GlobalSum; carefully in NORM_2 case
785  switch (type) {
786  case NORM_1_AND_2: {
787  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IceModelVec::norm_all(...): NORM_1_AND_2"
788  " not implemented (called as %s.norm_all(...))",
789  m_impl->name.c_str());
790 
791  }
792  case NORM_1: {
793  for (unsigned int k = 0; k < m_impl->dof; ++k) {
794  result[k] = GlobalSum(m_impl->grid->com, result[k]);
795  }
796  return result;
797  }
798  case NORM_2: {
799  for (unsigned int k = 0; k < m_impl->dof; ++k) {
800  // undo sqrt in VecNorm before sum; sum up; take sqrt
801  result[k] = sqrt(GlobalSum(m_impl->grid->com, result[k]*result[k]));
802  }
803  return result;
804  }
805  case NORM_INFINITY: {
806  for (unsigned int k = 0; k < m_impl->dof; ++k) {
807  result[k] = GlobalMax(m_impl->grid->com, result[k]);
808  }
809  return result;
810  }
811  default: {
812  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IceModelVec::norm_all(...): unknown norm type"
813  " (called as %s.norm_all(...))",
814  m_impl->name.c_str());
815  }
816  } // switch
817  } else {
818  return result;
819  }
820 }
821 
822 void IceModelVec::write(const std::string &filename) const {
823  // We expect the file to be present and ready to write into.
824  File file(m_impl->grid->com,
825  filename,
826  string_to_backend(m_impl->grid->ctx()->config()->get_string("output.format")),
828  m_impl->grid->ctx()->pio_iosys_id());
829 
830  this->write(file);
831 }
832 
833 void IceModelVec::read(const std::string &filename, unsigned int time) {
834  File file(m_impl->grid->com, filename, PISM_GUESS, PISM_READONLY);
835  this->read(file, time);
836 }
837 
838 void IceModelVec::regrid(const std::string &filename, RegriddingFlag flag,
839  double default_value) {
840  File file(m_impl->grid->com, filename, PISM_GUESS, PISM_READONLY);
841 
842  try {
843  this->regrid(file, flag, default_value);
844  } catch (RuntimeError &e) {
845  e.add_context("regridding '%s' from '%s'",
846  this->get_name().c_str(), filename.c_str());
847  throw;
848  }
849 }
850 
851 /** Read a field from a file, interpolating onto the current grid.
852  *
853  * When `flag` is set to `CRITICAL`, stop if could not find the variable
854  * in the provided input file; `default_value` is ignored.
855  *
856  * When `flag` is set to `OPTIONAL`, fill this IceModelVec with
857  * `default_value` if could not find the variable in the provided
858  * input file.
859  *
860  * When `flag` is set to `CRITICAL_FILL_MISSING`, replace missing
861  * values matching the `_FillValue` attribute with `default_value`,
862  * stop if could not find the variable.
863  *
864  * When `flag` is set to `OPTIONAL_FILL_MISSING`, replace missing
865  * values matching the `_FillValue` attribute with `default_value`;
866  * fill the whole IceModelVec with `default_value` if could not find
867  * the variable.
868  *
869  * @param file input file
870  * @param flag regridding mode, see above
871  * @param default_value default value, meaning depends on the
872  * regridding mode flag
873  *
874  * @return 0 on success
875  */
876 void IceModelVec::regrid(const File &file, RegriddingFlag flag,
877  double default_value) {
878  m_impl->grid->ctx()->log()->message(3, " [%s] Regridding %s...\n",
879  timestamp(m_impl->grid->com).c_str(), m_impl->name.c_str());
880  double start_time = get_time();
881  m_impl->grid->ctx()->profiling().begin("io.regridding");
882  {
883  this->regrid_impl(file, flag, default_value);
884  inc_state_counter(); // mark as modified
885  }
886  m_impl->grid->ctx()->profiling().end("io.regridding");
887  double
888  end_time = get_time(),
889  time_spent = end_time - start_time;
890 
891  if (time_spent > 1.0) {
892  m_impl->grid->ctx()->log()->message(3, " done in %f seconds.\n", time_spent);
893  } else {
894  m_impl->grid->ctx()->log()->message(3, " done.\n");
895  }
896 }
897 
898 void IceModelVec::read(const File &file, const unsigned int time) {
899  this->read_impl(file, time);
900  inc_state_counter(); // mark as modified
901 }
902 
903 void IceModelVec::write(const File &file) const {
904  define(file);
905 
906  m_impl->grid->ctx()->log()->message(3, " [%s] Writing %s...",
907  timestamp(m_impl->grid->com).c_str(),
908  m_impl->name.c_str());
909 
910  double start_time = get_time();
911  write_impl(file);
912  double end_time = get_time();
913 
914  const double
915  time_spent = end_time - start_time,
916  megabyte = pow(2, 20),
917  mb_double = sizeof(double) * size() / megabyte,
918  mb_float = sizeof(float) * size() / megabyte;
919 
920  std::string timestamp = pism::timestamp(m_impl->grid->com);
921  std::string spacer(timestamp.size(), ' ');
922  if (time_spent > 1) {
923  m_impl->grid->ctx()->log()->message(3,
924  "\n"
925  " [%s] Done writing %s (%f Mb double, %f Mb float)\n"
926  " %s in %f seconds (%f minutes).\n"
927  " %s Effective throughput: double: %f Mb/s, float: %f Mb/s.\n",
928  timestamp.c_str(), m_impl->name.c_str(), mb_double, mb_float,
929  spacer.c_str(), time_spent, time_spent / 60.0,
930  spacer.c_str(),
931  mb_double / time_spent, mb_float / time_spent);
932  } else {
933  m_impl->grid->ctx()->log()->message(3, " done.\n");
934  }
935 }
936 
938  // empty
939 }
940 
942  while (not m_vecs.empty()) {
943  try {
944  m_vecs.back()->end_access();
945  m_vecs.pop_back();
946  } catch (...) {
947  handle_fatal_errors(MPI_COMM_SELF);
948  }
949  }
950 }
951 
952 AccessList::AccessList(std::initializer_list<const PetscAccessible *> vecs) {
953  for (const auto *j : vecs) {
954  assert(j != nullptr);
955  add(*j);
956  }
957 }
958 
960  add(vec);
961 }
962 
964  vec.begin_access();
965  m_vecs.push_back(&vec);
966 }
967 
968 void AccessList::add(const std::vector<const PetscAccessible*> &vecs) {
969  for (const auto *v : vecs) {
970  assert(v != nullptr);
971  add(*v);
972  }
973 }
974 
975 //! Return the total number of elements in the *owned* part of an array.
976 size_t IceModelVec::size() const {
977  // m_impl->dof > 1 for vector, staggered grid 2D fields, etc. In this case
978  // zlevels.size() == 1. For 3D fields, m_impl->dof == 1 (all 3D fields are
979  // scalar) and zlevels.size() corresponds to dof of the underlying PETSc
980  // DM object.
981 
982  size_t
983  Mx = m_impl->grid->Mx(),
984  My = m_impl->grid->My(),
985  Mz = m_impl->zlevels.size(),
986  dof = m_impl->dof;
987 
988  return Mx * My * Mz * dof;
989 }
990 
991 /*! Allocate a copy on processor zero and the scatter needed to move data.
992  */
993 std::shared_ptr<petsc::Vec> IceModelVec::allocate_proc0_copy() const {
994  PetscErrorCode ierr;
995  Vec v_proc0 = NULL;
996  Vec result = NULL;
997 
998  ierr = PetscObjectQuery((PetscObject)dm()->get(), "v_proc0", (PetscObject*)&v_proc0);
999  PISM_CHK(ierr, "PetscObjectQuery")
1000  ;
1001  if (v_proc0 == NULL) {
1002 
1003  // natural_work will be destroyed at the end of scope, but it will
1004  // only decrement the reference counter incremented by
1005  // PetscObjectCompose below.
1006  petsc::Vec natural_work;
1007  // create a work vector with natural ordering:
1008  ierr = DMDACreateNaturalVector(*dm(), natural_work.rawptr());
1009  PISM_CHK(ierr, "DMDACreateNaturalVector");
1010 
1011  // this increments the reference counter of natural_work
1012  ierr = PetscObjectCompose((PetscObject)dm()->get(), "natural_work",
1013  (PetscObject)((::Vec)natural_work));
1014  PISM_CHK(ierr, "PetscObjectCompose");
1015 
1016  // scatter_to_zero will be destroyed at the end of scope, but it
1017  // will only decrement the reference counter incremented by
1018  // PetscObjectCompose below.
1019  petsc::VecScatter scatter_to_zero;
1020 
1021  // initialize the scatter to processor 0 and create storage on processor 0
1022  ierr = VecScatterCreateToZero(natural_work, scatter_to_zero.rawptr(),
1023  &v_proc0);
1024  PISM_CHK(ierr, "VecScatterCreateToZero");
1025 
1026  // this increments the reference counter of scatter_to_zero
1027  ierr = PetscObjectCompose((PetscObject)dm()->get(), "scatter_to_zero",
1028  (PetscObject)((::VecScatter)scatter_to_zero));
1029  PISM_CHK(ierr, "PetscObjectCompose");
1030 
1031  // this increments the reference counter of v_proc0
1032  ierr = PetscObjectCompose((PetscObject)dm()->get(), "v_proc0",
1033  (PetscObject)v_proc0);
1034  PISM_CHK(ierr, "PetscObjectCompose");
1035 
1036  // We DO NOT call VecDestroy(v_proc0): the petsc::Vec wrapper will
1037  // take care of this.
1038  result = v_proc0;
1039  } else {
1040  // We DO NOT call VecDestroy(result): the petsc::Vec wrapper will
1041  // take care of this.
1042  ierr = VecDuplicate(v_proc0, &result);
1043  PISM_CHK(ierr, "VecDuplicate");
1044  }
1045  return std::shared_ptr<petsc::Vec>(new petsc::Vec(result));
1046 }
1047 
1048 void IceModelVec::put_on_proc0(petsc::Vec &parallel, petsc::Vec &onp0) const {
1049  PetscErrorCode ierr = 0;
1050  VecScatter scatter_to_zero = NULL;
1051  Vec natural_work = NULL;
1052 
1053  ierr = PetscObjectQuery((PetscObject)dm()->get(), "scatter_to_zero",
1054  (PetscObject*)&scatter_to_zero);
1055  PISM_CHK(ierr, "PetscObjectQuery");
1056 
1057  ierr = PetscObjectQuery((PetscObject)dm()->get(), "natural_work",
1058  (PetscObject*)&natural_work);
1059  PISM_CHK(ierr, "PetscObjectQuery");
1060 
1061  if (natural_work == NULL || scatter_to_zero == NULL) {
1062  throw RuntimeError(PISM_ERROR_LOCATION, "call allocate_proc0_copy() before calling put_on_proc0");
1063  }
1064 
1065  ierr = DMDAGlobalToNaturalBegin(*dm(), parallel, INSERT_VALUES, natural_work);
1066  PISM_CHK(ierr, "DMDAGlobalToNaturalBegin");
1067 
1068  ierr = DMDAGlobalToNaturalEnd(*dm(), parallel, INSERT_VALUES, natural_work);
1069  PISM_CHK(ierr, "DMDAGlobalToNaturalEnd");
1070 
1071  ierr = VecScatterBegin(scatter_to_zero, natural_work, onp0,
1072  INSERT_VALUES, SCATTER_FORWARD);
1073  PISM_CHK(ierr, "VecScatterBegin");
1074 
1075  ierr = VecScatterEnd(scatter_to_zero, natural_work, onp0,
1076  INSERT_VALUES, SCATTER_FORWARD);
1077  PISM_CHK(ierr, "VecScatterEnd");
1078 }
1079 
1080 
1081 //! Puts a local IceModelVec2S on processor 0.
1083  if (m_impl->ghosted) {
1085  this->copy_to_vec(dm(), tmp);
1086  put_on_proc0(tmp, onp0);
1087  } else {
1088  put_on_proc0(vec(), onp0);
1089  }
1090 }
1091 
1092 void IceModelVec::get_from_proc0(petsc::Vec &onp0, petsc::Vec &parallel) const {
1093  PetscErrorCode ierr;
1094 
1095  VecScatter scatter_to_zero = NULL;
1096  Vec natural_work = NULL;
1097  ierr = PetscObjectQuery((PetscObject)dm()->get(), "scatter_to_zero",
1098  (PetscObject*)&scatter_to_zero);
1099  PISM_CHK(ierr, "PetscObjectQuery");
1100 
1101  ierr = PetscObjectQuery((PetscObject)dm()->get(), "natural_work",
1102  (PetscObject*)&natural_work);
1103  PISM_CHK(ierr, "PetscObjectQuery");
1104 
1105  if (natural_work == NULL || scatter_to_zero == NULL) {
1106  throw RuntimeError(PISM_ERROR_LOCATION, "call allocate_proc0_copy() before calling get_from_proc0");
1107  }
1108 
1109  ierr = VecScatterBegin(scatter_to_zero, onp0, natural_work,
1110  INSERT_VALUES, SCATTER_REVERSE);
1111  PISM_CHK(ierr, "VecScatterBegin");
1112 
1113  ierr = VecScatterEnd(scatter_to_zero, onp0, natural_work,
1114  INSERT_VALUES, SCATTER_REVERSE);
1115  PISM_CHK(ierr, "VecScatterEnd");
1116 
1117  ierr = DMDANaturalToGlobalBegin(*dm(), natural_work, INSERT_VALUES, parallel);
1118  PISM_CHK(ierr, "DMDANaturalToGlobalBegin");
1119 
1120  ierr = DMDANaturalToGlobalEnd(*dm(), natural_work, INSERT_VALUES, parallel);
1121  PISM_CHK(ierr, "DMDANaturalToGlobalEnd");
1122 }
1123 
1124 //! Gets a local IceModelVec2 from processor 0.
1126  if (m_impl->ghosted) {
1128  get_from_proc0(onp0, tmp);
1129  global_to_local(*dm(), tmp, vec());
1130  } else {
1131  get_from_proc0(onp0, vec());
1132  }
1133  inc_state_counter(); // mark as modified
1134 }
1135 
1136 /*!
1137  * Copy data to rank 0 and compute the checksum.
1138  *
1139  * Does not use ghosts. Results should be independent of the parallel domain
1140  * decomposition.
1141  */
1143 
1144  auto v = allocate_proc0_copy();
1145  put_on_proc0(*v);
1146 
1147  MPI_Comm com = m_impl->grid->ctx()->com();
1148 
1149  int rank = 0;
1150  MPI_Comm_rank(com, &rank);
1151 
1152  uint64_t result = 0;
1153  if (rank == 0) {
1154  petsc::VecArray array(*v);
1155 
1156  PetscInt size = 0;
1157  PetscErrorCode ierr = VecGetLocalSize(*v, &size);
1158  PISM_CHK(ierr, "VecGetLocalSize");
1159 
1160  result = pism::fletcher64((uint32_t*)array.get(), size * 2);
1161  }
1162  MPI_Bcast(&result, 1, MPI_UINT64_T, 0, com);
1163 
1164  return result;
1165 }
1166 
1167 /*!
1168  * Compute a checksum of a vector.
1169  *
1170  * The result depends on the number of processors used.
1171  *
1172  * We assume that sizeof(double) == 2 * sizeof(uint32_t), i.e. double uses 64 bits.
1173  */
1174 uint64_t IceModelVec::fletcher64() const {
1175  static_assert(sizeof(double) == 2 * sizeof(uint32_t), "Cannot compile IceModelVec::fletcher64() (sizeof(double) != 2 * sizeof(uint32_t))");
1176 
1177  MPI_Status mpi_stat;
1178  const int checksum_tag = 42;
1179 
1180  MPI_Comm com = m_impl->grid->ctx()->com();
1181 
1182  int rank = 0;
1183  MPI_Comm_rank(com, &rank);
1184 
1185  int comm_size = 0;
1186  MPI_Comm_size(com, &comm_size);
1187 
1188  PetscInt local_size = 0;
1189  PetscErrorCode ierr = VecGetLocalSize(vec(), &local_size); PISM_CHK(ierr, "VecGetLocalSize");
1190  uint64_t sum = 0;
1191  {
1192  petsc::VecArray v(vec());
1193  // compute checksums for local patches on all ranks
1194  sum = pism::fletcher64((uint32_t*)v.get(), local_size * 2);
1195  }
1196 
1197  if (rank == 0) {
1198  std::vector<uint64_t> sums(comm_size);
1199 
1200  // gather checksums of patches on rank 0
1201  sums[0] = sum;
1202  for (int r = 1; r < comm_size; ++r) {
1203  MPI_Recv(&sums[r], 1, MPI_UINT64_T, r, checksum_tag, com, &mpi_stat);
1204  }
1205 
1206  // compute the checksum of checksums
1207  sum = pism::fletcher64((uint32_t*)sums.data(), comm_size * 2);
1208  } else {
1209  MPI_Send(&sum, 1, MPI_UINT64_T, 0, checksum_tag, com);
1210  }
1211 
1212  // broadcast to all ranks
1213  MPI_Bcast(&sum, 1, MPI_UINT64_T, 0, com);
1214 
1215  return sum;
1216 }
1217 
1218 std::string IceModelVec::checksum(bool serial) const {
1219  if (serial) {
1220  // unsigned long long is supposed to be at least 64 bit long
1221  return pism::printf("%016llx", (unsigned long long int)this->fletcher64_serial());
1222  }
1223  // unsigned long long is supposed to be at least 64 bit long
1224  return pism::printf("%016llx", (unsigned long long int)this->fletcher64());
1225 }
1226 
1227 void IceModelVec::print_checksum(const char *prefix, bool serial) const {
1228  auto log = m_impl->grid->ctx()->log();
1229 
1230  log->message(1, "%s%s: %s\n", prefix, m_impl->name.c_str(), checksum(serial).c_str());
1231 }
1232 
1234  const std::string &spec1, const std::string &spec2) {
1235  units::Converter c(system, spec1, spec2);
1236 
1237  // has to be a PetscInt because of the VecGetLocalSize() call
1238  PetscInt data_size = 0;
1239  PetscErrorCode ierr = VecGetLocalSize(v, &data_size);
1240  PISM_CHK(ierr, "VecGetLocalSize");
1241 
1242  petsc::VecArray data(v);
1243  c.convert_doubles(data.get(), data_size);
1244 }
1245 
1247  const IceModelVec2Stag &input,
1248  bool include_floating_ice,
1249  IceModelVec2S &result) {
1250 
1251  using mask::grounded_ice;
1252  using mask::icy;
1253 
1254  assert(cell_type.stencil_width() > 0);
1255  assert(input.stencil_width() > 0);
1256 
1257  IceGrid::ConstPtr grid = result.grid();
1258 
1259  IceModelVec::AccessList list{&cell_type, &input, &result};
1260 
1261  for (Points p(*grid); p; p.next()) {
1262  const int i = p.i(), j = p.j();
1263 
1264  if (cell_type.grounded_ice(i, j) or
1265  (include_floating_ice and cell_type.icy(i, j))) {
1266  auto M = cell_type.star(i, j);
1267  auto F = input.star(i, j);
1268 
1269  int n = 0, e = 0, s = 0, w = 0;
1270  if (include_floating_ice) {
1271  n = static_cast<int>(icy(M.n));
1272  e = static_cast<int>(icy(M.e));
1273  s = static_cast<int>(icy(M.s));
1274  w = static_cast<int>(icy(M.w));
1275  } else {
1276  n = static_cast<int>(grounded_ice(M.n));
1277  e = static_cast<int>(grounded_ice(M.e));
1278  s = static_cast<int>(grounded_ice(M.s));
1279  w = static_cast<int>(grounded_ice(M.w));
1280  }
1281 
1282  if (n + e + s + w > 0) {
1283  result(i, j) = (n * F.n + e * F.e + s * F.s + w * F.w) / (n + e + s + w);
1284  } else {
1285  result(i, j) = 0.0;
1286  }
1287  } else {
1288  result(i, j) = 0.0;
1289  }
1290  }
1291 }
1292 
1294  const IceModelVec2Stag &input,
1295  bool include_floating_ice,
1296  IceModelVec2V &result) {
1297 
1298  using mask::grounded_ice;
1299  using mask::icy;
1300 
1301  assert(cell_type.stencil_width() > 0);
1302  assert(input.stencil_width() > 0);
1303 
1304  IceGrid::ConstPtr grid = result.grid();
1305 
1306  IceModelVec::AccessList list{&cell_type, &input, &result};
1307 
1308  for (Points p(*grid); p; p.next()) {
1309  const int i = p.i(), j = p.j();
1310 
1311  auto M = cell_type.star(i, j);
1312  auto F = input.star(i, j);
1313 
1314  int n = 0, e = 0, s = 0, w = 0;
1315  if (include_floating_ice) {
1316  n = static_cast<int>(icy(M.n));
1317  e = static_cast<int>(icy(M.e));
1318  s = static_cast<int>(icy(M.s));
1319  w = static_cast<int>(icy(M.w));
1320  } else {
1321  n = static_cast<int>(grounded_ice(M.n));
1322  e = static_cast<int>(grounded_ice(M.e));
1323  s = static_cast<int>(grounded_ice(M.s));
1324  w = static_cast<int>(grounded_ice(M.w));
1325  }
1326 
1327  if (e + w > 0) {
1328  result(i, j).u = (e * F.e + w * F.w) / (e + w);
1329  } else {
1330  result(i, j).u = 0.0;
1331  }
1332 
1333  if (n + s > 0) {
1334  result(i, j).v = (n * F.n + s * F.s) / (n + s);
1335  } else {
1336  result(i, j).v = 0.0;
1337  }
1338  }
1339 }
1340 
1341 //! \brief View a 2D vector field using existing PETSc viewers.
1342 void IceModelVec::view(std::vector<std::shared_ptr<petsc::Viewer> > viewers) const {
1343  PetscErrorCode ierr;
1344 
1345  if (ndims() == 3) {
1347  "cannot 'view' a 3D field '%s'",
1348  get_name().c_str());
1349  }
1350 
1351  // Get the dof=1, stencil_width=0 DMDA (components are always scalar
1352  // and we just need a global Vec):
1353  auto da2 = m_impl->grid->get_dm(1, 0);
1354 
1355  petsc::TemporaryGlobalVec tmp(da2);
1356 
1357  for (unsigned int i = 0; i < ndof(); ++i) {
1358  std::string
1359  long_name = m_impl->metadata[i].get_string("long_name"),
1360  units = m_impl->metadata[i].get_string("units"),
1361  glaciological_units = m_impl->metadata[i].get_string("glaciological_units"),
1362  title = long_name + " (" + glaciological_units + ")";
1363 
1364  PetscViewer v = *viewers[i].get();
1365 
1366  PetscDraw draw = NULL;
1367  ierr = PetscViewerDrawGetDraw(v, 0, &draw);
1368  PISM_CHK(ierr, "PetscViewerDrawGetDraw");
1369 
1370  ierr = PetscDrawSetTitle(draw, title.c_str());
1371  PISM_CHK(ierr, "PetscDrawSetTitle");
1372 
1373  get_dof(da2, tmp, i);
1374 
1375  convert_vec(tmp, m_impl->metadata[i].unit_system(),
1376  units, glaciological_units);
1377 
1378  double bounds[2] = {0.0, 0.0};
1379  ierr = VecMin(tmp, NULL, &bounds[0]); PISM_CHK(ierr, "VecMin");
1380  ierr = VecMax(tmp, NULL, &bounds[1]); PISM_CHK(ierr, "VecMax");
1381 
1382  if (bounds[0] == bounds[1]) {
1383  bounds[0] = -1.0;
1384  bounds[1] = 1.0;
1385  }
1386 
1387  ierr = PetscViewerDrawSetBounds(v, 1, bounds);
1388  PISM_CHK(ierr, "PetscViewerDrawSetBounds");
1389 
1390  ierr = VecView(tmp, v);
1391  PISM_CHK(ierr, "VecView");
1392  }
1393 }
1394 
1395 } // end of namespace pism
std::vector< const PetscAccessible * > m_vecs
Definition: iceModelVec.hh:68
void add(const PetscAccessible &v)
Definition: iceModelVec.cc:963
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool icy(int i, int j) const
bool grounded_ice(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
stencils::Star< int > star(int i, int j) const
stencils::Star< double > star(int i, int j) const
Returns the values at interfaces of the cell i,j using the staggered grid.
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
std::vector< int > shape() const
Definition: iceModelVec.cc:156
size_t size() const
Return the total number of elements in the owned part of an array.
Definition: iceModelVec.cc:976
void regrid_impl(const File &file, RegriddingFlag flag, double default_value=0.0)
Gets an IceModelVec from a file file, interpolating onto the current grid.
Definition: iceModelVec.cc:422
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
void get_dof(std::shared_ptr< petsc::DM > da_result, petsc::Vec &result, unsigned int start, unsigned int count=1) const
Definition: iceModelVec.cc:277
void read_impl(const File &file, unsigned int time)
Reads appropriate NetCDF variable(s) into an IceModelVec.
Definition: iceModelVec.cc:476
std::string checksum(bool serial) const
void set_name(const std::string &name)
Sets the variable name to name.
Definition: iceModelVec.cc:374
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:334
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void view(std::vector< std::shared_ptr< petsc::Viewer > > viewers) const
View a 2D vector field using existing PETSc viewers.
void dump(const char filename[]) const
Dumps a variable to a file, overwriting this file's contents (for debugging).
Definition: iceModelVec.cc:583
uint64_t fletcher64() const
uint64_t fletcher64_serial() const
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
Definition: iceModelVec.cc:838
void print_checksum(const char *prefix="", bool serial=false) const
virtual void begin_access() const
Checks if an IceModelVec is allocated and calls DAVecGetArray.
Definition: iceModelVec.cc:622
petsc::Vec & vec() const
Definition: iceModelVec.cc:342
std::shared_ptr< petsc::DM > dm() const
Definition: iceModelVec.cc:356
virtual ~IceModelVec()
Definition: iceModelVec.cc:101
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void get_from_proc0(petsc::Vec &onp0)
Gets a local IceModelVec2 from processor 0.
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
void read(const std::string &filename, unsigned int time)
Definition: iceModelVec.cc:833
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void put_on_proc0(petsc::Vec &onp0) const
Puts a local IceModelVec2S on processor 0.
void check_array_indices(int i, int j, unsigned int k) const
Check array indices and warn if they are out of range.
Definition: iceModelVec.cc:690
void set_dof(std::shared_ptr< petsc::DM > da_source, petsc::Vec &source, unsigned int start, unsigned int count=1)
Definition: iceModelVec.cc:304
std::array< double, 2 > range() const
Result: min <- min(v[j]), max <- max(v[j]).
Definition: iceModelVec.cc:193
std::vector< double > levels() const
Definition: iceModelVec.cc:136
void checkCompatibility(const char *function, const IceModelVec &other) const
Checks if two IceModelVecs have compatible sizes, dimensions and numbers of degrees of freedom.
Definition: iceModelVec.cc:599
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
Definition: iceModelVec.cc:523
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
virtual void end_access() const
Checks if an IceModelVec is allocated and calls DAVecRestoreArray.
Definition: iceModelVec.cc:643
void copy_to_vec(std::shared_ptr< petsc::DM > destination_da, petsc::Vec &destination) const
Copies v to a global vector 'destination'. Ghost points are discarded.
Definition: iceModelVec.cc:266
void write_impl(const File &file) const
Writes an IceModelVec to a NetCDF file.
Definition: iceModelVec.cc:544
std::shared_ptr< petsc::Vec > allocate_proc0_copy() const
Definition: iceModelVec.cc:993
void set_begin_access_use_dof(bool flag)
Definition: iceModelVec.cc:181
std::vector< double > norm(int n) const
Computes the norm of all the components of an IceModelVec.
Definition: iceModelVec.cc:769
unsigned int ndims() const
Returns the number of spatial dimensions.
Definition: iceModelVec.cc:152
const std::string & get_name() const
Get the name of an IceModelVec object.
Definition: iceModelVec.cc:386
void inc_state_counter()
Increment the object state counter.
Definition: iceModelVec.cc:147
int state_counter() const
Get the object state counter.
Definition: iceModelVec.cc:124
void shift(double alpha)
Result: v[j] <- v[j] + alpha for all j. Calls VecShift.
Definition: iceModelVec.cc:244
unsigned int ndof() const
Returns the number of degrees of freedom per grid point.
Definition: iceModelVec.cc:132
IceModelVec(IceGrid::ConstPtr grid, const std::string &name, IceModelVecKind ghostedp, size_t dof, size_t stencil_width, const std::vector< double > &zlevels)
Definition: iceModelVec.cc:56
void add(double alpha, const IceModelVec &x)
Result: v <- v + alpha * x. Calls VecAXPY.
Definition: iceModelVec.cc:234
void set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
Definition: iceModelVec.cc:175
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Definition: iceModelVec.hh:202
std::shared_ptr< const Logger > ConstPtr
Definition: Logger.hh:46
void failed()
Indicates a failure of a parallel section.
virtual void begin_access() const =0
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
IO_Type get_output_type() const
double * get()
Definition: Vec.cc:54
Wrapper around VecGetArray and VecRestoreArray.
Definition: Vec.hh:44
void convert_doubles(double *data, size_t length) const
Definition: Units.cc:252
std::shared_ptr< System > Ptr
Definition: Units.hh:47
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
void append_time(const File &file, const Config &config, double time_seconds)
Prepare a file for output.
Definition: io_helpers.cc:256
void read_spatial_variable(const SpatialVariableMetadata &variable, const IceGrid &grid, const File &file, unsigned int time, double *output)
Read a variable from a file into an array output.
Definition: io_helpers.cc:616
void define_spatial_variable(const SpatialVariableMetadata &var, const IceGrid &grid, const File &file, IO_Type default_type)
Define a NetCDF variable corresponding to a VariableMetadata object.
Definition: io_helpers.cc:558
void write_spatial_variable(const SpatialVariableMetadata &var, const IceGrid &grid, const File &file, const double *input)
Write a double array to a file.
Definition: io_helpers.cc:719
void regrid_spatial_variable(SpatialVariableMetadata &var, const IceGrid &grid, const File &file, RegriddingFlag flag, bool report_range, bool allow_extrapolation, double default_value, InterpolationType interpolation_type, double *output)
Regrid from a NetCDF file into a distributed array output.
Definition: io_helpers.cc:783
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
Definition: io_helpers.cc:219
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:47
bool grounded_ice(int M)
Definition: Mask.hh:50
void compute_params(const IceModelVec *const x, const IceModelVec *const y, const IceModelVec *const z, int &ghosts, bool &scatter)
Compute parameters for 2D loop computations involving 3 IceModelVecs.
Definition: iceModelVec.cc:728
double get_time()
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static double F(double SL, double B, double H, double alpha)
IO_Type
Definition: IO_Flags.hh:31
@ PISM_NAT
Definition: IO_Flags.hh:32
@ PISM_DOUBLE
Definition: IO_Flags.hh:38
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 GlobalMax(MPI_Comm comm, double *local, double *result, int count)
IO_Backend string_to_backend(const std::string &backend)
Definition: File.cc:60
static NormType int_to_normtype(int input)
Definition: iceModelVec.cc:219
std::string printf(const char *format,...)
uint64_t fletcher64(const uint32_t *data, size_t length)
static const double k
Definition: exactTestP.cc:45
static double end_time(const Config &config, double time_start, const std::string &calendar, const units::Unit &time_units)
Definition: Time.cc:398
RegriddingFlag
Definition: IO_Flags.hh:70
static void global_to_local(petsc::DM &dm, Vec source, Vec destination)
Definition: iceModelVec.cc:46
@ PISM_GUESS
Definition: IO_Flags.hh:41
std::string timestamp(MPI_Comm com)
Creates a time-stamp used for the history NetCDF attribute.
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
void handle_fatal_errors(MPI_Comm com)
@ PISM_READWRITE_CLOBBER
create a file for writing, overwrite if present
Definition: IO_Flags.hh:53
@ PISM_READWRITE
open an existing file for reading and writing
Definition: IO_Flags.hh:51
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
void GlobalMin(MPI_Comm comm, double *local, double *result, int count)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
void convert_vec(petsc::Vec &v, units::System::Ptr system, const std::string &spec1, const std::string &spec2)
void staggered_to_regular(const IceModelVec2CellType &cell_type, const IceModelVec2Stag &input, bool include_floating_ice, IceModelVec2S &result)
IceModelVecKind
What "kind" of a vector to create: with or without ghosts.
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
double sum(const IceModelVec2S &input)
Sums up all the values in an IceModelVec2S object. Ignores ghosts.
bool report_range
If true, report range when regridding.
InterpolationType interpolation_type
std::shared_ptr< petsc::DM > da
distributed mesh manager (DM)
std::vector< double > zlevels
Vertical levels (for 3D fields)
unsigned int da_stencil_width
stencil width supported by the DA
int state_counter
Internal IceModelVec "revision number".
std::vector< SpatialVariableMetadata > metadata
Metadata (NetCDF variable attributes)
bool begin_access_use_dof
If true, use DMDAVecGetArrayDOF() in begin_access()
unsigned int dof
number of "degrees of freedom" per grid point
gsl_interp_accel * bsearch_accel
bool ghosted
true if this IceModelVec is ghosted
IceGrid::ConstPtr grid
The computational grid.
int count
Definition: test_cube.c:16