PISM, A Parallel Ice Sheet Model  stable v1.2 committed by Constantine Khrulev on 2020-02-11 20:24:05 -0900
PicoGeometry.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2019 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 <algorithm> // max_element
21 
22 #include "PicoGeometry.hh"
23 #include "pism/util/connected_components.hh"
24 #include "pism/util/IceModelVec2CellType.hh"
25 #include "pism/util/pism_utilities.hh"
26 
27 namespace pism {
28 namespace ocean {
29 
31  : Component(grid),
32  m_continental_shelf(grid, "pico_contshelf_mask", WITHOUT_GHOSTS),
33  m_boxes(grid, "pico_box_mask", WITHOUT_GHOSTS),
34  m_ice_shelves(grid, "pico_shelf_mask", WITHOUT_GHOSTS),
35  m_distance_gl(grid, "pico_distance_gl", WITH_GHOSTS),
36  m_distance_cf(grid, "pico_distance_cf", WITH_GHOSTS),
37  m_ocean_mask(grid, "pico_ocean_mask", WITH_GHOSTS),
38  m_lake_mask(grid, "pico_lake_mask", WITHOUT_GHOSTS),
39  m_ice_rises(grid, "pico_ice_rise_mask", WITH_GHOSTS),
40  m_tmp(grid, "temporary_storage", WITHOUT_GHOSTS) {
41 
42  m_boxes.metadata().set_number("_FillValue", 0.0);
43 
44  m_ice_rises.metadata().set_numbers("flag_values",
46  m_ice_rises.metadata().set_string("flag_meanings",
47  "ocean ice_rise continental_ice_sheet, floating_ice");
48 
50 }
51 
53  // empty
54 }
55 
57  return m_continental_shelf;
58 }
59 
61  return m_boxes;
62 }
63 
65  return m_ice_shelves;
66 }
67 
69  return m_ice_rises;
70 }
71 
72 /*!
73  * Compute masks needed by the PICO physics code.
74  *
75  * After this call box_mask(), ice_shelf_mask(), and continental_shelf_mask() will be up
76  * to date.
77  */
78 void PicoGeometry::update(const IceModelVec2S &bed_elevation, const IceModelVec2CellType &cell_type) {
79  bool exclude_ice_rises = m_config->get_flag("ocean.pico.exclude_ice_rises");
80 
81  int n_boxes = m_config->get_number("ocean.pico.number_of_boxes");
82 
83  double continental_shelf_depth = m_config->get_number("ocean.pico.continental_shelf_depth");
84 
85  // these three could be done at the same time
86  {
87  compute_ice_rises(cell_type, exclude_ice_rises, m_ice_rises);
88 
89  compute_ocean_mask(cell_type, m_ocean_mask);
90 
91  compute_lakes(cell_type, m_lake_mask);
92  }
93 
94  // these two could be optimized by trying to reduce the number of times we update ghosts
95  {
98 
100 
102  }
103 
104  // these two could be done at the same time
105  {
107 
108  compute_continental_shelf_mask(bed_elevation, m_ice_rises, continental_shelf_depth, m_continental_shelf);
109  }
110 
112 }
113 
114 
116 
117 /*!
118  * Re-label components in a mask processed by label_connected_components.
119  *
120  * If type is `BY_AREA`, the biggest one gets the value of 2, all the other ones 1, the
121  * background is set to zero.
122  *
123  * If type is `AREA_THRESHOLD`, patches with areas above `threshold` get the value of 2,
124  * all the other ones 1, the background is set to zero.
125  */
126 static void relabel(RelabelingType type,
127  double threshold,
128  IceModelVec2Int &mask) {
129 
130  IceGrid::ConstPtr grid = mask.grid();
131 
132  int max_index = mask.range().max;
133 
134  if (max_index < 1) {
135  // No components labeled. Fill the mask with zeros and quit.
136  mask.set(0.0);
137  return;
138  }
139 
140  std::vector<double> area(max_index + 1, 0.0);
141  {
142 
143  ParallelSection loop(grid->com);
144  try {
145  for (Points p(*grid); p; p.next()) {
146  const int i = p.i(), j = p.j();
147 
148  int index = mask(i, j);
149 
150  if (index > max_index or index < 0) {
151  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid component index: %d", index);
152  }
153 
154  if (index > 0) {
155  // count areas of actual components, ignoring the background (index == 0)
156  area[index] += 1.0;
157  }
158  }
159  } catch (...) {
160  loop.failed();
161  }
162  loop.check();
163 
164  for (unsigned int k = 0; k < area.size(); ++k) {
165  area[k] = grid->cell_area() * GlobalSum(grid->com, area[k]);
166  }
167  }
168 
169  if (type == BY_AREA) {
170  // find the biggest component
171  int biggest_component = 0;
172  for (unsigned int k = 0; k < area.size(); ++k) {
173  if (area[k] > area[biggest_component]) {
174  biggest_component = k;
175  }
176  }
177 
178  // re-label
179  for (Points p(*grid); p; p.next()) {
180  const int i = p.i(), j = p.j();
181 
182  int component_index = mask.as_int(i, j);
183 
184  if (component_index == biggest_component) {
185  mask(i, j) = 2.0;
186  } else if (component_index > 0) {
187  mask(i, j) = 1.0;
188  } else {
189  mask(i, j) = 0.0;
190  }
191  }
192  } else {
193  for (Points p(*grid); p; p.next()) {
194  const int i = p.i(), j = p.j();
195 
196  int component_index = mask.as_int(i, j);
197 
198  if (area[component_index] > threshold) {
199  mask(i, j) = 2.0;
200  } else if (component_index > 0) {
201  mask(i, j) = 1.0;
202  } else {
203  mask(i, j) = 0.0;
204  }
205  }
206  }
207 }
208 
209 /*!
210  * Run the serial connected-component labeling algorithm on m_tmp.
211  */
214 
215  ParallelSection rank0(m_grid->com);
216  try {
217  if (m_grid->rank() == 0) {
218  petsc::VecArray mask_p0(*m_tmp_p0);
219  label_connected_components(mask_p0.get(), m_grid->My(), m_grid->Mx(), false, 0.0);
220  }
221  } catch (...) {
222  rank0.failed();
223  }
224  rank0.check();
225 
227 }
228 
229 static bool edge_p(int i, int j, int Mx, int My) {
230  return (i == 0) or (i == Mx - 1) or (j == 0) or (j == My - 1);
231 }
232 
233 /*!
234  * Compute the mask identifying "subglacial lakes", i.e. floating ice areas that are not
235  * connected to the open ocean.
236  *
237  * Resulting mask contains:
238  *
239  * 0 - grounded ice
240  * 1 - floating ice not connected to the open ocean
241  * 2 - floating ice or ice-free ocean connected to the open ocean
242  */
244  IceModelVec::AccessList list{ &cell_type, &m_tmp };
245 
246  const int
247  Mx = m_grid->Mx(),
248  My = m_grid->My();
249 
250  // assume that ocean points (i.e. floating, either icy or ice-free) at the edge of the
251  // domain belong to the "open ocean"
252  for (Points p(*m_grid); p; p.next()) {
253  const int i = p.i(), j = p.j();
254 
255  if (cell_type.ocean(i, j)) {
256  m_tmp(i, j) = 1.0;
257 
258  if (edge_p(i, j, Mx, My)) {
259  m_tmp(i, j) = 2.0;
260  }
261  } else {
262  m_tmp(i, j) = 0.0;
263  }
264  }
265 
266  // identify "floating" areas that are not connected to the open ocean as defined above
267  {
269 
270  ParallelSection rank0(m_grid->com);
271  try {
272  if (m_grid->rank() == 0) {
273  petsc::VecArray mask_p0(*m_tmp_p0);
274  label_connected_components(mask_p0.get(), My, Mx, true, 2.0);
275  }
276  } catch (...) {
277  rank0.failed();
278  }
279  rank0.check();
280 
282  }
283 
284  result.copy_from(m_tmp);
285 }
286 
287 /*!
288  * Compute the mask identifying ice rises, i.e. grounded ice areas not connected to the
289  * continental ice sheet.
290  *
291  * Resulting mask contains:
292  *
293  * 0 - ocean
294  * 1 - ice rises
295  * 2 - continental ice sheet
296  * 3 - floating ice
297  */
298 void PicoGeometry::compute_ice_rises(const IceModelVec2CellType &cell_type, bool exclude_ice_rises,
299  IceModelVec2Int &result) {
300  IceModelVec::AccessList list{ &cell_type, &m_tmp };
301 
302  // mask of zeros and ones: one if grounded ice, zero otherwise
303  for (Points p(*m_grid); p; p.next()) {
304  const int i = p.i(), j = p.j();
305 
306  if (cell_type.grounded(i, j)) {
307  m_tmp(i, j) = 2.0;
308  } else {
309  m_tmp(i, j) = 0.0;
310  }
311  }
312 
313  if (exclude_ice_rises) {
314  label_tmp();
315 
317  m_config->get_number("ocean.pico.maximum_ice_rise_area", "m2"),
318  m_tmp);
319  }
320 
321  // mark floating ice areas in this mask (reduces the number of masks we need later)
322  for (Points p(*m_grid); p; p.next()) {
323  const int i = p.i(), j = p.j();
324 
325  if (m_tmp(i, j) == 0.0 and cell_type.icy(i, j)) {
326  m_tmp(i, j) = FLOATING;
327  }
328  }
329 
330  result.copy_from(m_tmp);
331 }
332 
333 /*!
334  * Compute the continental ice shelf mask.
335  *
336  * Resulting mask contains:
337  *
338  * 0 - ocean or icy
339  * 1 - ice-free areas with bed elevation > threshold and not connected to the continental ice sheet
340  * 2 - ice-free areas with bed elevation > threshold, connected to the continental ice sheet
341  */
343  const IceModelVec2Int &ice_rise_mask, double bed_elevation_threshold,
344  IceModelVec2Int &result) {
345  IceModelVec::AccessList list{ &bed_elevation, &ice_rise_mask, &m_tmp };
346 
347  for (Points p(*m_grid); p; p.next()) {
348  const int i = p.i(), j = p.j();
349 
350  m_tmp(i, j) = 0.0;
351 
352  if (bed_elevation(i, j) > bed_elevation_threshold) {
353  m_tmp(i, j) = 1.0;
354  }
355 
356  if (ice_rise_mask.as_int(i, j) == CONTINENTAL) {
357  m_tmp(i, j) = 2.0;
358  }
359  }
360 
361  // use "iceberg identification" to label parts *not* connected to the continental ice
362  // sheet
363  {
365 
366  ParallelSection rank0(m_grid->com);
367  try {
368  if (m_grid->rank() == 0) {
369  petsc::VecArray mask_p0(*m_tmp_p0);
370  label_connected_components(mask_p0.get(), m_grid->My(), m_grid->Mx(), true, 2.0);
371  }
372  } catch (...) {
373  rank0.failed();
374  }
375  rank0.check();
376 
378  }
379 
380  // At this point areas with bed > threshold are 1, everything else is zero.
381  //
382  // Now we need to mark the continental shelf itself.
383  for (Points p(*m_grid); p; p.next()) {
384  const int i = p.i(), j = p.j();
385 
386  if (m_tmp(i, j) > 0.0) {
387  continue;
388  }
389 
390  if (bed_elevation(i, j) > bed_elevation_threshold and
391  ice_rise_mask.as_int(i, j) == OCEAN) {
392  m_tmp(i, j) = 2.0;
393  }
394  }
395 
396  result.copy_from(m_tmp);
397 }
398 
399 /*!
400  * Compute the mask identifying ice shelves.
401  *
402  * Each shelf gets an individual integer label.
403  *
404  * Two shelves connected by an ice rise are considered to be parts of the same shelf.
405  *
406  * Floating ice cells that are not connected to the ocean ("subglacial lakes") are
407  * excluded.
408  */
409 void PicoGeometry::compute_ice_shelf_mask(const IceModelVec2Int &ice_rise_mask, const IceModelVec2Int &lake_mask,
410  IceModelVec2Int &result) {
411  IceModelVec::AccessList list{ &ice_rise_mask, &lake_mask, &m_tmp };
412 
413  for (Points p(*m_grid); p; p.next()) {
414  const int i = p.i(), j = p.j();
415 
416  int M = ice_rise_mask.as_int(i, j);
417 
418  if (M == RISE or M == FLOATING) {
419  m_tmp(i, j) = 1.0;
420  } else {
421  m_tmp(i, j) = 0.0;
422  }
423  }
424 
425  label_tmp();
426 
427  // remove ice rises and lakes
428  for (Points p(*m_grid); p; p.next()) {
429  const int i = p.i(), j = p.j();
430 
431  if (ice_rise_mask.as_int(i, j) == RISE or lake_mask.as_int(i, j) == 1) {
432  m_tmp(i, j) = 0.0;
433  }
434  }
435 
436  result.copy_from(m_tmp);
437 }
438 
439 /*!
440  * Compute the mask identifying ice-free ocean and "holes" in ice shelves.
441  *
442  * Resulting mask contains:
443  *
444  * - 0 - icy cells
445  * - 1 - ice-free ocean which is not connected to the open ocean
446  * - 2 - open ocean
447  *
448  */
450  IceModelVec::AccessList list{ &cell_type, &m_tmp };
451 
452  // mask of zeros and ones: one if ice-free ocean, zero otherwise
453  for (Points p(*m_grid); p; p.next()) {
454  const int i = p.i(), j = p.j();
455 
456  if (cell_type.ice_free_ocean(i, j)) {
457  m_tmp(i, j) = 1.0;
458  } else {
459  m_tmp(i, j) = 0.0;
460  }
461  }
462 
463  label_tmp();
464 
465  relabel(BY_AREA, 0.0, m_tmp);
466 
467  result.copy_from(m_tmp);
468 }
469 
471  const IceModelVec2Int &ice_rises,
472  bool exclude_ice_rises,
473  IceModelVec2Int &result) {
474 
475  IceModelVec::AccessList list{ &ice_rises, &ocean_mask, &result };
476 
477  result.set(-1);
478 
479  // Find the grounding line and the ice front and set result to 1 if ice shelf cell is
480  // next to the grounding line, Ice holes within the shelf are treated like ice shelf
481  // cells, if exclude_ice_rises is set then ice rises are also treated as ice shelf cells.
482 
483  ParallelSection loop(m_grid->com);
484  try {
485  for (Points p(*m_grid); p; p.next()) {
486  const int i = p.i(), j = p.j();
487 
488  if (ice_rises.as_int(i, j) == FLOATING or
489  ocean_mask.as_int(i, j) == 1 or
490  (exclude_ice_rises and ice_rises.as_int(i, j) == RISE)) {
491  // if this is an ice shelf cell (or an ice rise) or a hole in an ice shelf
492 
493  // label the shelf cells adjacent to the grounding line with 1
494  bool neighbor_to_land =
495  (ice_rises(i, j + 1) == CONTINENTAL or
496  ice_rises(i, j - 1) == CONTINENTAL or
497  ice_rises(i + 1, j) == CONTINENTAL or
498  ice_rises(i - 1, j) == CONTINENTAL or
499  ice_rises(i + 1, j + 1) == CONTINENTAL or
500  ice_rises(i + 1, j - 1) == CONTINENTAL or
501  ice_rises(i - 1, j + 1) == CONTINENTAL or
502  ice_rises(i - 1, j - 1) == CONTINENTAL);
503 
504  if (neighbor_to_land) {
505  // i.e. there is a grounded neighboring cell (which is not an ice rise!)
506  result(i, j) = 1;
507  } else {
508  result(i, j) = 0;
509  }
510  }
511  }
512  } catch (...) {
513  loop.failed();
514  }
515  loop.check();
516 
517  result.update_ghosts();
518 
519  eikonal_equation(result);
520 }
521 
523  const IceModelVec2Int &ice_rises,
524  bool exclude_ice_rises,
525  IceModelVec2Int &result) {
526 
527  IceModelVec::AccessList list{ &ice_rises, &ocean_mask, &result };
528 
529  result.set(-1);
530 
531  ParallelSection loop(m_grid->com);
532  try {
533  for (Points p(*m_grid); p; p.next()) {
534  const int i = p.i(), j = p.j();
535 
536  if (ice_rises.as_int(i, j) == FLOATING or
537  ocean_mask.as_int(i, j) == 1 or
538  (exclude_ice_rises and ice_rises.as_int(i, j) == RISE)) {
539  // if this is an ice shelf cell (or an ice rise) or a hole in an ice shelf
540 
541  // label the shelf cells adjacent to the ice front with 1
542  auto M = ocean_mask.int_star(i, j);
543 
544  if (M.n == 2 or M.e == 2 or M.s == 2 or M.w == 2) {
545  // i.e. there is a neighboring open ocean cell
546  result(i, j) = 1;
547  } else {
548  result(i, j) = 0;
549  }
550  }
551  }
552  } catch (...) {
553  loop.failed();
554  }
555  loop.check();
556 
557  result.update_ghosts();
558 
559  eikonal_equation(result);
560 }
561 
562 /*!
563  * Find an approximate solution of the Eikonal equation on a given domain.
564  *
565  * To specify the problem, the input field (mask) should be filled with
566  *
567  * - negative values outside the domain
568  * - zeros within the domain
569  * - ones at "wave front" locations
570  *
571  * For example, to compute distances from the grounding line within ice shelves, fill
572  * generic ice shelf locations with zeros, set neighbors of the grounding line to 1, and
573  * the rest of the grid with -1 or some other negative number.
574  *
575  * Note that this implementation updates ghosts *every* iteration. We could speed this
576  * up by checking if a point at a boundary of the processor sub-domain was updated and
577  * update ghosts in those cases only.
578  */
580 
581  assert(mask.stencil_width() > 0);
582 
583  IceGrid::ConstPtr grid = mask.grid();
584 
585  double current_label = 1;
586  double continue_loop = 1;
587  while (continue_loop != 0) {
588 
589  continue_loop = 0;
590 
591  for (Points p(*grid); p; p.next()) {
592  const int i = p.i(), j = p.j();
593 
594  if (mask.as_int(i, j) == 0) {
595 
596  auto R = mask.int_star(i, j);
597 
598  if (R.ij == 0 and
599  (R.n == current_label or R.s == current_label or
600  R.e == current_label or R.w == current_label)) {
601  // i.e. this is an shelf cell with no distance assigned yet and with a neighbor
602  // that has a distance assigned
603  mask(i, j) = current_label + 1;
604  continue_loop = 1;
605  }
606  }
607  } // loop over grid points
608 
609  current_label++;
610  mask.update_ghosts();
611 
612  continue_loop = GlobalMax(grid->com, continue_loop);
613  }
614 }
615 
617  const IceModelVec2Int &shelf_mask, int max_number_of_boxes,
618  IceModelVec2Int &result) {
619 
620  IceModelVec::AccessList list{ &D_gl, &D_cf, &shelf_mask, &result };
621 
622  int n_shelves = shelf_mask.range().max + 1;
623 
624  std::vector<double> GL_distance_max(n_shelves, 0.0);
625  std::vector<double> CF_distance_max(n_shelves, 0.0);
626 
627  for (Points p(*m_grid); p; p.next()) {
628  const int i = p.i(), j = p.j();
629 
630  int shelf_id = shelf_mask(i, j);
631  assert(shelf_id >= 0);
632  assert(shelf_id < n_shelves + 1);
633 
634  if (shelf_id == 0) {
635  // not at a shelf; skip to the next grid point
636  continue;
637  }
638 
639  if (D_gl(i, j) > GL_distance_max[shelf_id]) {
640  GL_distance_max[shelf_id] = D_gl(i, j);
641  }
642 
643  if (D_cf(i, j) > CF_distance_max[shelf_id]) {
644  CF_distance_max[shelf_id] = D_cf(i, j);
645  }
646  }
647 
648  // compute global maximums
649  for (int k = 0; k < n_shelves; ++k) {
650  GL_distance_max[k] = GlobalMax(m_grid->com, GL_distance_max[k]);
651  CF_distance_max[k] = GlobalMax(m_grid->com, CF_distance_max[k]);
652  }
653 
654  double GL_distance_ref = *std::max_element(GL_distance_max.begin(), GL_distance_max.end());
655 
656  // compute the number of boxes in each shelf
657 
658  std::vector<int> n_boxes(n_shelves, 0);
659  int n_min = 1;
660  double zeta = 0.5;
661 
662  for (int k = 0; k < n_shelves; ++k) {
663  n_boxes[k] = n_min + round(pow((GL_distance_max[k] / GL_distance_ref), zeta) * (max_number_of_boxes - n_min));
664 
665  n_boxes[k] = std::min(n_boxes[k], max_number_of_boxes);
666  }
667 
668  result.set(0.0);
669 
670  for (Points p(*m_grid); p; p.next()) {
671  const int i = p.i(), j = p.j();
672 
673  int d_gl = D_gl.as_int(i, j);
674  int d_cf = D_cf.as_int(i, j);
675 
676  if (shelf_mask.as_int(i, j) > 0 and d_gl > 0.0 and d_cf > 0.0) {
677  int shelf_id = shelf_mask(i, j);
678  int n = n_boxes[shelf_id];
679 
680  // relative position on the shelf (ranges from 0 to 1), increasing towards the
681  // calving front (see equation 10 in the PICO paper)
682  double r = d_gl / (double)(d_gl + d_cf);
683 
684  double C = pow(1.0 - r, 2.0);
685  for (int k = 0; k < n; ++k) {
686  if ((n - k - 1) / (double)n <= C and C <= (n - k) / (double)n) {
687  result(i, j) = std::min(d_gl, k + 1);
688  }
689  }
690  }
691  }
692 }
693 
694 } // end of namespace ocean
695 } // end of namespace pism
pism::GlobalMax
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
Definition: pism_utilities.cc:138
n
#define n
Definition: exactTestM.c:37
pism::Range::max
double max
Definition: iceModelVec.hh:49
pism::ocean::PicoGeometry::compute_ocean_mask
void compute_ocean_mask(const IceModelVec2CellType &cell_type, IceModelVec2Int &result)
Definition: PicoGeometry.cc:449
pism::ocean::PicoGeometry::m_distance_cf
IceModelVec2Int m_distance_cf
Definition: PicoGeometry.hh:78
PISM_ERROR_LOCATION
#define PISM_ERROR_LOCATION
Definition: error_handling.hh:42
pism::IceModelVec::update_ghosts
virtual void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:630
pism::ocean::PicoGeometry::ice_rise_mask
const IceModelVec2Int & ice_rise_mask() const
Definition: PicoGeometry.cc:68
pism::Points
Definition: IceGrid.hh:391
pism::PointsWithGhosts::next
void next()
Definition: IceGrid.hh:363
pism::ocean::PicoGeometry::compute_ice_shelf_mask
void compute_ice_shelf_mask(const IceModelVec2Int &ice_rises_mask, const IceModelVec2Int &lake_mask, IceModelVec2Int &result)
Definition: PicoGeometry.cc:409
pism::ocean::BY_AREA
@ BY_AREA
Definition: PicoGeometry.cc:115
pism::IceModelVec::set
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:685
pism::IceGrid::ConstPtr
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:215
pism::ocean::PicoGeometry::m_ice_shelves
IceModelVec2Int m_ice_shelves
Definition: PicoGeometry.hh:74
pism::ocean::PicoGeometry::CONTINENTAL
@ CONTINENTAL
Definition: PicoGeometry.hh:50
pism
Definition: AgeColumnSystem.cc:24
pism::ocean::PicoGeometry::m_boxes
IceModelVec2Int m_boxes
Definition: PicoGeometry.hh:73
pism::Component::m_config
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:140
pism::RuntimeError::formatted
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Definition: error_handling.cc:47
pism::ocean::PicoGeometry::compute_distances_cf
void compute_distances_cf(const IceModelVec2Int &ocean_mask, const IceModelVec2Int &ice_rises, bool exclude_ice_rises, IceModelVec2Int &dist_cf)
Definition: PicoGeometry.cc:522
pism::IceModelVec2S::copy_from
virtual void copy_from(const IceModelVec &source)
Result: v <- source. Leaves metadata alone but copies values in Vec. Uses VecCopy.
Definition: iceModelVec2.cc:492
pism::Component::m_grid
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:138
pism::IceModelVec::allocate_proc0_copy
petsc::Vec::Ptr allocate_proc0_copy() const
Definition: iceModelVec.cc:1100
pism::IceModelVec2CellType::ocean
bool ocean(int i, int j) const
Definition: IceModelVec2CellType.hh:45
pism::VariableMetadata::set_string
void set_string(const std::string &name, const std::string &value)
Set a string attribute.
Definition: VariableMetadata.cc:343
pism::ocean::AREA_THRESHOLD
@ AREA_THRESHOLD
Definition: PicoGeometry.cc:115
pism::ocean::PicoGeometry::compute_distances_gl
void compute_distances_gl(const IceModelVec2Int &ocean_mask, const IceModelVec2Int &ice_rises, bool exclude_ice_rises, IceModelVec2Int &dist_gl)
Definition: PicoGeometry.cc:470
pism::ocean::PicoGeometry::m_ice_rises
IceModelVec2Int m_ice_rises
Definition: PicoGeometry.hh:81
pism::IceModelVec2CellType::ice_free_ocean
bool ice_free_ocean(int i, int j) const
Definition: IceModelVec2CellType.hh:69
pism::ocean::PicoGeometry::continental_shelf_mask
const IceModelVec2Int & continental_shelf_mask() const
Definition: PicoGeometry.cc:56
pism::IceModelVec::stencil_width
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:357
pism::ocean::PicoGeometry::box_mask
const IceModelVec2Int & box_mask() const
Definition: PicoGeometry.cc:60
pism::ocean::eikonal_equation
void eikonal_equation(IceModelVec2Int &mask)
Definition: PicoGeometry.cc:579
pism::k
static const double k
Definition: exactTestP.cc:45
double
pism::ocean::PicoGeometry::m_lake_mask
IceModelVec2Int m_lake_mask
Definition: PicoGeometry.hh:80
pism::ocean::PicoGeometry::update
void update(const IceModelVec2S &bed_elevation, const IceModelVec2CellType &cell_type)
Definition: PicoGeometry.cc:78
pism::AccessList
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:67
pism::ParallelSection
Definition: error_handling.hh:70
pism::VariableMetadata::set_number
void set_number(const std::string &name, double value)
Set a scalar attribute to a single (scalar) value.
Definition: VariableMetadata.cc:292
pism::ocean::PicoGeometry::m_continental_shelf
IceModelVec2Int m_continental_shelf
Definition: PicoGeometry.hh:72
pism::ocean::PicoGeometry::compute_box_mask
void compute_box_mask(const IceModelVec2Int &D_gl, const IceModelVec2Int &D_cf, const IceModelVec2Int &shelf_mask, int n_boxes, IceModelVec2Int &result)
Definition: PicoGeometry.cc:616
pism::IceModelVec2Int
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:488
pism::WITH_GHOSTS
@ WITH_GHOSTS
Definition: iceModelVec.hh:46
pism::ocean::PicoGeometry::m_tmp_p0
petsc::Vec::Ptr m_tmp_p0
Definition: PicoGeometry.hh:85
pism::ParallelSection::check
void check()
Definition: error_handling.cc:205
label_connected_components
void label_connected_components(double *image, unsigned int n_rows, unsigned int n_cols, bool identify_icebergs, double mask_grounded)
In-place labeling of connected components using a 2-scan algorithm with run-length encoding.
Definition: connected_components.cc:47
pism::IceModelVec2S
Definition: iceModelVec.hh:435
pism::ocean::PicoGeometry::PicoGeometry
PicoGeometry(IceGrid::ConstPtr grid)
Definition: PicoGeometry.cc:30
pism::ocean::PicoGeometry::label_tmp
void label_tmp()
Definition: PicoGeometry.cc:212
C
const int C[]
Definition: ssafd_code.cc:44
pism::mask::ocean
bool ocean(int M)
An ocean cell (floating ice or ice-free).
Definition: Mask.hh:39
pism::IceModelVec2CellType::grounded
bool grounded(int i, int j) const
Definition: IceModelVec2CellType.hh:49
pism::ocean::PicoGeometry::m_distance_gl
IceModelVec2Int m_distance_gl
Definition: PicoGeometry.hh:77
pism::IceModelVec2Int::int_star
StarStencil< int > int_star(int i, int j) const
Definition: IceModelVec_inline.hh:103
pism::ocean::PicoGeometry::RISE
@ RISE
Definition: PicoGeometry.hh:50
pism::ocean::PicoGeometry::FLOATING
@ FLOATING
Definition: PicoGeometry.hh:50
pism::IceModelVec::grid
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:79
pism::IceModelVec2CellType::icy
bool icy(int i, int j) const
Definition: IceModelVec2CellType.hh:53
pism::IceModelVec::put_on_proc0
void put_on_proc0(Vec onp0) const
Puts a local IceModelVec2S on processor 0.
Definition: iceModelVec.cc:1189
pism::petsc::VecArray::get
double * get()
Definition: Vec.cc:54
pism::ocean::edge_p
static bool edge_p(int i, int j, int Mx, int My)
Definition: PicoGeometry.cc:229
pism::WITHOUT_GHOSTS
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:46
pism::VariableMetadata::set_numbers
void set_numbers(const std::string &name, const std::vector< double > &values)
Set a scalar attribute to a single (scalar) value.
Definition: VariableMetadata.cc:297
pism::ocean::RelabelingType
RelabelingType
Definition: PicoGeometry.cc:115
pism::ocean::PicoGeometry::compute_ice_rises
void compute_ice_rises(const IceModelVec2CellType &cell_type, bool exclude_ice_rises, IceModelVec2Int &result)
Definition: PicoGeometry.cc:298
pism::ocean::PicoGeometry::ice_shelf_mask
const IceModelVec2Int & ice_shelf_mask() const
Definition: PicoGeometry.cc:64
pism::ocean::PicoGeometry::compute_lakes
void compute_lakes(const IceModelVec2CellType &cell_type, IceModelVec2Int &result)
Definition: PicoGeometry.cc:243
pism::ocean::PicoGeometry::OCEAN
@ OCEAN
Definition: PicoGeometry.hh:50
pism::IceModelVec2CellType
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
Definition: IceModelVec2CellType.hh:29
pism::ocean::PicoGeometry::m_ocean_mask
IceModelVec2Int m_ocean_mask
Definition: PicoGeometry.hh:79
pism::IceModelVec::get_from_proc0
void get_from_proc0(Vec onp0)
Gets a local IceModelVec2 from processor 0.
Definition: iceModelVec.cc:1232
pism::IceModelVec2Int::as_int
int as_int(int i, int j) const
Definition: IceModelVec_inline.hh:95
pism::ocean::PicoGeometry::compute_continental_shelf_mask
void compute_continental_shelf_mask(const IceModelVec2S &bed_elevation, const IceModelVec2Int &ice_rises_mask, double bed_elevation_threshold, IceModelVec2Int &result)
Definition: PicoGeometry.cc:342
pism::petsc::VecArray
Wrapper around VecGetArray and VecRestoreArray.
Definition: Vec.hh:44
pism::IceModelVec::range
virtual Range range() const
Result: min <- min(v[j]), max <- max(v[j]).
Definition: iceModelVec.cc:151
pism::ParallelSection::failed
void failed()
Indicates a failure of a parallel section.
Definition: error_handling.cc:189
pism::GlobalSum
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
Definition: pism_utilities.cc:142
pism::ocean::PicoGeometry::m_tmp
IceModelVec2Int m_tmp
Definition: PicoGeometry.hh:84
pism::IceModelVec::metadata
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:507
pism::ocean::PicoGeometry::~PicoGeometry
virtual ~PicoGeometry()
Definition: PicoGeometry.cc:52
pism::Component
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:103
pism::ocean::relabel
static void relabel(RelabelingType type, double threshold, IceModelVec2Int &mask)
Definition: PicoGeometry.cc:126
PicoGeometry.hh