PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
PicoGeometry.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2019, 2020, 2021 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 #include "PicoGeometry.hh"
22 #include "pism/util/connected_components.hh"
23 #include "pism/util/IceModelVec2CellType.hh"
24 #include "pism/util/pism_utilities.hh"
25 #include "pism/util/petscwrappers/Vec.hh"
26 
27 #include "pism/coupler/util/options.hh"
28 
29 namespace pism {
30 namespace ocean {
31 
33  : Component(grid),
34  m_continental_shelf(grid, "pico_contshelf_mask", WITHOUT_GHOSTS),
35  m_boxes(grid, "pico_box_mask", WITHOUT_GHOSTS),
36  m_ice_shelves(grid, "pico_shelf_mask", WITHOUT_GHOSTS),
37  m_basin_mask(m_grid, "basins", WITH_GHOSTS),
38  m_distance_gl(grid, "pico_distance_gl", WITH_GHOSTS),
39  m_distance_cf(grid, "pico_distance_cf", WITH_GHOSTS),
40  m_ocean_mask(grid, "pico_ocean_mask", WITH_GHOSTS),
41  m_lake_mask(grid, "pico_lake_mask", WITHOUT_GHOSTS),
42  m_ice_rises(grid, "pico_ice_rise_mask", WITH_GHOSTS),
43  m_tmp(grid, "temporary_storage", WITHOUT_GHOSTS) {
44 
45  m_boxes.metadata()["_FillValue"] = {0.0};
46 
47  m_ice_rises.metadata()["flag_values"] = {OCEAN, RISE, CONTINENTAL, FLOATING};
48  m_ice_rises.metadata()["flag_meanings"] =
49  "ocean ice_rise continental_ice_sheet, floating_ice";
50 
51  m_basin_mask.set_attrs("climate_forcing", "mask determines basins for PICO",
52  "", "", "", 0);
53  m_n_basins = 0;
54 
56 }
57 
59  return m_continental_shelf;
60 }
61 
63  return m_boxes;
64 }
65 
67  return m_ice_shelves;
68 }
69 
71  return m_ice_rises;
72 }
73 
75  return m_basin_mask;
76 }
77 
79 
80  ForcingOptions opt(*m_grid->ctx(), "ocean.pico");
81 
83 
84  m_n_basins = static_cast<int>(max(m_basin_mask)) + 1;
85 }
86 
87 /*!
88  * Compute masks needed by the PICO physics code.
89  *
90  * After this call box_mask(), ice_shelf_mask(), and continental_shelf_mask() will be up
91  * to date.
92  */
93 void PicoGeometry::update(const IceModelVec2S &bed_elevation,
94  const IceModelVec2CellType &cell_type) {
95 
96  // Update basin adjacency.
97  //
98  // basin_neighbors() below uses the cell type mask to find
99  // adjacent basins by iterating over the current ice front. This means that basin
100  // adjacency cannot be pre-computed during initialization.
101  {
103 
104  // report
105  for (const auto &p : m_basin_neighbors) {
106  std::vector<std::string> neighbors;
107  for (const auto &n : p.second) {
108  neighbors.emplace_back(pism::printf("%d", n));
109  }
110  std::string neighbor_list = pism::join(neighbors, ", ");
111  m_log->message(3, "PICO: basin %d neighbors: %s\n",
112  p.first, neighbor_list.c_str());
113  }
114  }
115 
116  bool exclude_ice_rises = m_config->get_flag("ocean.pico.exclude_ice_rises");
117 
118  // these three could be done at the same time
119  {
120  compute_ice_rises(cell_type, exclude_ice_rises, m_ice_rises);
121 
122  compute_ocean_mask(cell_type, m_ocean_mask);
123 
124  compute_lakes(cell_type, m_lake_mask);
125 
126  }
127 
128  // these two could be optimized by trying to reduce the number of times we update ghosts
129  {
132 
134 
136  }
137 
138  // computing ice_shelf_mask and box_mask could be done at the same time
139  {
141  auto n_shelves = static_cast<int>(max(m_ice_shelves)) + 1;
142 
143  std::vector<int> cfs_in_basins_per_shelf(n_shelves*m_n_basins, 0);
144  std::vector<int> most_shelf_cells_in_basin(n_shelves, 0);
146  most_shelf_cells_in_basin, cfs_in_basins_per_shelf);
147 
149  most_shelf_cells_in_basin, cfs_in_basins_per_shelf, n_shelves,
150  m_ice_shelves);
151 
152  double continental_shelf_depth = m_config->get_number("ocean.pico.continental_shelf_depth");
153 
154  compute_continental_shelf_mask(bed_elevation, m_ice_rises, continental_shelf_depth,
156  }
157 
158  int n_boxes = static_cast<int>(m_config->get_number("ocean.pico.number_of_boxes"));
159 
161 }
162 
163 
165 
166 /*!
167  * Re-label components in a mask processed by label_connected_components.
168  *
169  * If type is `BY_AREA`, the biggest one gets the value of 2, all the other ones 1, the
170  * background is set to zero.
171  *
172  * If type is `AREA_THRESHOLD`, patches with areas above `threshold` get the value of 2,
173  * all the other ones 1, the background is set to zero.
174  */
175 static void relabel(RelabelingType type,
176  double threshold,
177  IceModelVec2Int &mask) {
178 
179  IceGrid::ConstPtr grid = mask.grid();
180 
181  int max_index = static_cast<int>(mask.range()[1]);
182 
183  if (max_index < 1) {
184  // No components labeled. Fill the mask with zeros and quit.
185  mask.set(0.0);
186  return;
187  }
188 
189  std::vector<double> area(max_index + 1, 0.0);
190  std::vector<double> area1(max_index + 1, 0.0);
191  {
192 
193  ParallelSection loop(grid->com);
194  try {
195  for (Points p(*grid); p; p.next()) {
196  const int i = p.i(), j = p.j();
197 
198  int index = mask.as_int(i, j);
199 
200  if (index > max_index or index < 0) {
201  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid component index: %d", index);
202  }
203 
204  if (index > 0) {
205  // count areas of actual components, ignoring the background (index == 0)
206  area[index] += 1.0;
207  }
208  }
209  } catch (...) {
210  loop.failed();
211  }
212  loop.check();
213  GlobalSum(grid->com, area.data(), area1.data(), area.size());
214 
215  // copy data
216  area = area1;
217 
218  for (unsigned int k = 0; k < area.size(); ++k) {
219  area[k] = grid->cell_area() * area[k];
220  }
221  }
222 
223  if (type == BY_AREA) {
224  // find the biggest component
225  int biggest_component = 0;
226  for (unsigned int k = 0; k < area.size(); ++k) {
227  if (area[k] > area[biggest_component]) {
228  biggest_component = static_cast<int>(k);
229  }
230  }
231 
232  // re-label
233  for (Points p(*grid); p; p.next()) {
234  const int i = p.i(), j = p.j();
235 
236  int component_index = mask.as_int(i, j);
237 
238  if (component_index == biggest_component) {
239  mask(i, j) = 2.0;
240  } else if (component_index > 0) {
241  mask(i, j) = 1.0;
242  } else {
243  mask(i, j) = 0.0;
244  }
245  }
246  } else {
247  for (Points p(*grid); p; p.next()) {
248  const int i = p.i(), j = p.j();
249 
250  int component_index = mask.as_int(i, j);
251 
252  if (area[component_index] > threshold) {
253  mask(i, j) = 2.0;
254  } else if (component_index > 0) {
255  mask(i, j) = 1.0;
256  } else {
257  mask(i, j) = 0.0;
258  }
259  }
260  }
261 }
262 
263 /*!
264  * Run the serial connected-component labeling algorithm on m_tmp.
265  */
268 
269  ParallelSection rank0(m_grid->com);
270  try {
271  if (m_grid->rank() == 0) {
272  petsc::VecArray mask_p0(*m_tmp_p0);
274  static_cast<int>(m_grid->My()),
275  static_cast<int>(m_grid->Mx()),
276  false,
277  0.0);
278  }
279  } catch (...) {
280  rank0.failed();
281  }
282  rank0.check();
283 
285 }
286 
287 /*!
288  * Compute the mask identifying "subglacial lakes", i.e. floating ice areas that are not
289  * connected to the open ocean.
290  *
291  * Resulting mask contains:
292  *
293  * 0 - grounded ice
294  * 1 - floating ice not connected to the open ocean
295  * 2 - floating ice or ice-free ocean connected to the open ocean
296  */
298  IceModelVec::AccessList list{ &cell_type, &m_tmp };
299 
300  const int
301  Mx = m_grid->Mx(),
302  My = m_grid->My();
303 
304  // assume that ocean points (i.e. floating, either icy or ice-free) at the edge of the
305  // domain belong to the "open ocean"
306  for (Points p(*m_grid); p; p.next()) {
307  const int i = p.i(), j = p.j();
308 
309  if (cell_type.ocean(i, j)) {
310  m_tmp(i, j) = 1.0;
311 
312  if (grid_edge(*m_grid, i, j)) {
313  m_tmp(i, j) = 2.0;
314  }
315  } else {
316  m_tmp(i, j) = 0.0;
317  }
318  }
319 
320  // identify "floating" areas that are not connected to the open ocean as defined above
321  {
323 
324  ParallelSection rank0(m_grid->com);
325  try {
326  if (m_grid->rank() == 0) {
327  petsc::VecArray mask_p0(*m_tmp_p0);
328  label_connected_components(mask_p0.get(), My, Mx, true, 2.0);
329  }
330  } catch (...) {
331  rank0.failed();
332  }
333  rank0.check();
334 
336  }
337 
338  result.copy_from(m_tmp);
339 }
340 
341 /*!
342  * Compute the mask identifying ice rises, i.e. grounded ice areas not connected to the
343  * continental ice sheet.
344  *
345  * Resulting mask contains:
346  *
347  * 0 - ocean
348  * 1 - ice rises
349  * 2 - continental ice sheet
350  * 3 - floating ice
351  */
352 void PicoGeometry::compute_ice_rises(const IceModelVec2CellType &cell_type, bool exclude_ice_rises,
353  IceModelVec2Int &result) {
354  IceModelVec::AccessList list{ &cell_type, &m_tmp };
355 
356  // mask of zeros and ones: one if grounded ice, zero otherwise
357  for (Points p(*m_grid); p; p.next()) {
358  const int i = p.i(), j = p.j();
359 
360  if (cell_type.grounded(i, j)) {
361  m_tmp(i, j) = 2.0;
362  } else {
363  m_tmp(i, j) = 0.0;
364  }
365  }
366 
367  if (exclude_ice_rises) {
368  label_tmp();
369 
371  m_config->get_number("ocean.pico.maximum_ice_rise_area", "m2"),
372  m_tmp);
373  }
374 
375  // mark floating ice areas in this mask (reduces the number of masks we need later)
376  for (Points p(*m_grid); p; p.next()) {
377  const int i = p.i(), j = p.j();
378 
379  if (m_tmp(i, j) == 0.0 and cell_type.icy(i, j)) {
380  m_tmp(i, j) = FLOATING;
381  }
382  }
383 
384  result.copy_from(m_tmp);
385 }
386 
387 /*!
388  * Compute the continental ice shelf mask.
389  *
390  * Resulting mask contains:
391  *
392  * 0 - ocean or icy
393  * 1 - ice-free areas with bed elevation > threshold and not connected to the continental ice sheet
394  * 2 - ice-free areas with bed elevation > threshold, connected to the continental ice sheet
395  */
397  const IceModelVec2Int &ice_rise_mask,
398  double bed_elevation_threshold,
399  IceModelVec2Int &result) {
400  IceModelVec::AccessList list{ &bed_elevation, &ice_rise_mask, &m_tmp };
401 
402  for (Points p(*m_grid); p; p.next()) {
403  const int i = p.i(), j = p.j();
404 
405  m_tmp(i, j) = 0.0;
406 
407  if (bed_elevation(i, j) > bed_elevation_threshold) {
408  m_tmp(i, j) = 1.0;
409  }
410 
411  if (ice_rise_mask.as_int(i, j) == CONTINENTAL) {
412  m_tmp(i, j) = 2.0;
413  }
414  }
415 
416  // use "iceberg identification" to label parts *not* connected to the continental ice
417  // sheet
418  {
420 
421  ParallelSection rank0(m_grid->com);
422  try {
423  if (m_grid->rank() == 0) {
424  petsc::VecArray mask_p0(*m_tmp_p0);
425  label_connected_components(mask_p0.get(), m_grid->My(), m_grid->Mx(), true, 2.0);
426  }
427  } catch (...) {
428  rank0.failed();
429  }
430  rank0.check();
431 
433  }
434 
435  // At this point areas with bed > threshold are 1, everything else is zero.
436  //
437  // Now we need to mark the continental shelf itself.
438  for (Points p(*m_grid); p; p.next()) {
439  const int i = p.i(), j = p.j();
440 
441  if (m_tmp(i, j) > 0.0) {
442  continue;
443  }
444 
445  if (bed_elevation(i, j) > bed_elevation_threshold and
446  ice_rise_mask.as_int(i, j) == OCEAN) {
447  m_tmp(i, j) = 2.0;
448  }
449  }
450 
451  result.copy_from(m_tmp);
452 }
453 
454 /*!
455  * Compute the mask identifying ice shelves.
456  *
457  * Each shelf gets an individual integer label.
458  *
459  * Two shelves connected by an ice rise are considered to be parts of the same shelf.
460  *
461  * Floating ice cells that are not connected to the ocean ("subglacial lakes") are
462  * excluded.
463  */
464 void PicoGeometry::compute_ice_shelf_mask(const IceModelVec2Int &ice_rise_mask, const IceModelVec2Int &lake_mask,
465  IceModelVec2Int &result) {
466  IceModelVec::AccessList list{ &ice_rise_mask, &lake_mask, &m_tmp };
467 
468  for (Points p(*m_grid); p; p.next()) {
469  const int i = p.i(), j = p.j();
470 
471  int M = ice_rise_mask.as_int(i, j);
472 
473  if (M == RISE or M == FLOATING) {
474  m_tmp(i, j) = 1.0;
475  } else {
476  m_tmp(i, j) = 0.0;
477  }
478  }
479 
480  label_tmp();
481 
482  // remove ice rises and lakes
483  for (Points p(*m_grid); p; p.next()) {
484  const int i = p.i(), j = p.j();
485 
486  if (ice_rise_mask.as_int(i, j) == RISE or lake_mask.as_int(i, j) == 1) {
487  m_tmp(i, j) = 0.0;
488  }
489  }
490 
491  result.copy_from(m_tmp);
492 }
493 
494 /*!
495  * Compute the mask identifying ice-free ocean and "holes" in ice shelves.
496  *
497  * Resulting mask contains:
498  *
499  * - 0 - icy cells
500  * - 1 - ice-free ocean which is not connected to the open ocean
501  * - 2 - open ocean
502  *
503  */
505  IceModelVec::AccessList list{ &cell_type, &m_tmp };
506 
507  // mask of zeros and ones: one if ice-free ocean, zero otherwise
508  for (Points p(*m_grid); p; p.next()) {
509  const int i = p.i(), j = p.j();
510 
511  if (cell_type.ice_free_ocean(i, j)) {
512  m_tmp(i, j) = 1.0;
513  } else {
514  m_tmp(i, j) = 0.0;
515  }
516  }
517 
518  label_tmp();
519 
520  relabel(BY_AREA, 0.0, m_tmp);
521 
522  result.copy_from(m_tmp);
523 }
524 
525 /*!
526  * Find neighboring basins by checking for the basin boundaries on the ice free ocean.
527  *
528  * Should we identify the intersection at the coastline instead?
529  *
530  * Returns the map from the basin index to a set of indexes of neighbors.
531  */
532 std::map<int,std::set<int> > PicoGeometry::basin_neighbors(const IceModelVec2CellType &cell_type,
533  const IceModelVec2Int &basin_mask) {
534  using mask::ice_free_ocean;
535 
536  // Allocate the adjacency matrix. This uses twice the amount of storage necessary (the
537  // matrix is symmetric), but in known cases (i.e. with around 20 basins) we're wasting
538  // only ~200*sizeof(int) bytes, which is not that bad (and the code is a bit simpler).
539  std::vector<int> adjacency_matrix(m_n_basins * m_n_basins, 0);
540 
541  // short-cuts
542  auto mark_as_neighbors = [&](int b1, int b2) {
543  adjacency_matrix[b1 * m_n_basins + b2] = 1;
544  // preserve symmetry:
545  adjacency_matrix[b2 * m_n_basins + b1] = 1;
546  };
547 
548  auto adjacent = [&](int b1, int b2) {
549  return adjacency_matrix[b1 * m_n_basins + b2] > 0;
550  };
551 
552  IceModelVec::AccessList list{ &cell_type, &basin_mask };
553 
554  for (Points p(*m_grid); p; p.next()) {
555  const int i = p.i(), j = p.j();
556 
557  auto B = basin_mask.star(i, j);
558 
559  bool next_to_icefront = (cell_type.ice_free_ocean(i, j) and cell_type.next_to_ice(i,j));
560 
561  // skip the "dummy" basin and cells that are not at the ice front
562  if (B.ij == 0 or not next_to_icefront) {
563  continue;
564  }
565 
566  // Zero out IDs of basins for cell neighbors that are outside the modeling domain.
567  //
568  // This prevents "wrap-around" at grid boundaries.
569  {
570  B.n *= static_cast<int>(j < (int)m_grid->My() - 1);
571  B.e *= static_cast<int>(i < (int)m_grid->Mx() - 1);
572  B.s *= static_cast<int>(j > 0);
573  B.w *= static_cast<int>(i > 0);
574  }
575 
576  auto M = cell_type.star(i, j);
577 
578  if (ice_free_ocean(M.n)) {
579  mark_as_neighbors(B.ij, B.n);
580  }
581 
582  if (ice_free_ocean(M.s)) {
583  mark_as_neighbors(B.ij, B.s);
584  }
585 
586  if (ice_free_ocean(M.e)) {
587  mark_as_neighbors(B.ij, B.e);
588  }
589 
590  if (ice_free_ocean(M.w)) {
591  mark_as_neighbors(B.ij, B.w);
592  }
593  }
594 
595  // Make a copy to allreduce efficiently with IntelMPI:
596  {
597  std::vector<int> tmp(adjacency_matrix.size(), 0);
598  GlobalMax(m_grid->com, adjacency_matrix.data(), tmp.data(),
599  static_cast<int>(tmp.size()));
600  // Copy results:
601  adjacency_matrix = tmp;
602  }
603 
604  // Convert the matrix into a map "basin ID -> set of neighbors' IDs":
605  std::map<int,std::set<int> > result;
606  for (int b1 = 1; b1 < m_n_basins; ++b1) {
607  for (int b2 = b1 + 1; b2 < m_n_basins; ++b2) {
608  if (adjacent(b1, b2)) {
609  result[b1].insert(b2);
610  result[b2].insert(b1);
611  }
612  }
613  }
614 
615  return result;
616 }
617 
618 /*!
619  * Find all basins b, in which the ice shelf s has a calving front with potential ocean water intrusion.
620  * Find the basin bmax, in which the ice shelf s has the most cells
621  */
623  const IceModelVec2Int &basin_mask,
624  const IceModelVec2Int &shelf_mask,
625  int n_shelves,
626  std::vector<int> &most_shelf_cells_in_basin,
627  std::vector<int> &cfs_in_basins_per_shelf) {
628 
629  std::vector<int> n_shelf_cells_per_basin(n_shelves * m_n_basins,0);
630  // additional vectors to allreduce efficiently with IntelMPI
631  std::vector<int> n_shelf_cells_per_basinr(n_shelves * m_n_basins,0);
632  std::vector<int> cfs_in_basins_per_shelfr(n_shelves * m_n_basins,0);
633  std::vector<int> most_shelf_cells_in_basinr(n_shelves, 0);
634 
635  IceModelVec::AccessList list{ &cell_type, &basin_mask, &shelf_mask };
636 
637  {
638  for (Points p(*m_grid); p; p.next()) {
639  const int i = p.i(), j = p.j();
640  int s = shelf_mask.as_int(i, j);
641  int b = basin_mask.as_int(i, j);
642  int sb = s * m_n_basins + b;
643  n_shelf_cells_per_basin[sb]++;
644 
645  if (cell_type.as_int(i, j) == MASK_FLOATING) {
646  auto M = cell_type.star(i, j);
647  if (M.n == MASK_ICE_FREE_OCEAN or
648  M.e == MASK_ICE_FREE_OCEAN or
649  M.s == MASK_ICE_FREE_OCEAN or
650  M.w == MASK_ICE_FREE_OCEAN) {
651  if (cfs_in_basins_per_shelf[sb] != b) {
652  cfs_in_basins_per_shelf[sb] = b;
653  }
654  }
655  }
656  }
657 
658  GlobalSum(m_grid->com, cfs_in_basins_per_shelf.data(),
659  cfs_in_basins_per_shelfr.data(), n_shelves*m_n_basins);
660  GlobalSum(m_grid->com, n_shelf_cells_per_basin.data(),
661  n_shelf_cells_per_basinr.data(), n_shelves*m_n_basins);
662  // copy values
663  cfs_in_basins_per_shelf = cfs_in_basins_per_shelfr;
664  n_shelf_cells_per_basin = n_shelf_cells_per_basinr;
665 
666  for (int s = 0; s < n_shelves; s++) {
667  int n_shelf_cells_per_basin_max = 0;
668  for (int b = 0; b < m_n_basins; b++) {
669  int sb = s * m_n_basins + b;
670  if (n_shelf_cells_per_basin[sb] > n_shelf_cells_per_basin_max) {
671  most_shelf_cells_in_basin[s] = b;
672  n_shelf_cells_per_basin_max = n_shelf_cells_per_basin[sb];
673  }
674  }
675  }
676  }
677 }
678 
679 
680 /*!
681  * Find all ice shelves s that spread across non-neighboring basins with calving fronts in those basins and add an ice shelf mask number.
682  */
684  const IceModelVec2Int &basin_mask,
685  const std::map<int, std::set<int> > &basin_neighbors,
686  const std::vector<int> &most_shelf_cells_in_basin,
687  const std::vector<int> &cfs_in_basins_per_shelf,
688  int n_shelves,
689  IceModelVec2Int &shelf_mask) {
690  m_tmp.copy_from(shelf_mask);
691 
692  std::vector<int> n_shelf_cells_to_split(n_shelves * m_n_basins, 0);
693 
694  IceModelVec::AccessList list{ &cell_type, &basin_mask, &shelf_mask, &m_tmp };
695 
696  ParallelSection loop(m_grid->com);
697  try {
698  for (Points p(*m_grid); p; p.next()) {
699  const int i = p.i(), j = p.j();
700  if (cell_type.as_int(i, j) == MASK_FLOATING) {
701  int b = basin_mask.as_int(i, j);
702  int s = shelf_mask.as_int(i, j);
703  int b0 = most_shelf_cells_in_basin[s];
704  // basin_neighbors.at(b) may throw
705  bool neighbors = basin_neighbors.at(b).count(b0) > 0;
706  if (b != b0 and (not neighbors) and
707  cfs_in_basins_per_shelf[s * m_n_basins + b] > 0) {
708  n_shelf_cells_to_split[s * m_n_basins + b]++;
709  }
710  }
711  }
712  } catch (...) {
713  loop.failed();
714  }
715  loop.check();
716 
717  {
718  std::vector<int> tmp(n_shelves * m_n_basins, 0);
719  GlobalSum(m_grid->com, n_shelf_cells_to_split.data(),
720  tmp.data(), n_shelves * m_n_basins);
721  // copy values
722  n_shelf_cells_to_split = tmp;
723  }
724 
725  // no GlobalSum needed here, only local:
726  std::vector<int> add_shelf_instance(n_shelves * m_n_basins, 0);
727  int n_shelf_numbers_to_add = 0;
728  for (int s = 0; s < n_shelves; s++) {
729  int b0 = most_shelf_cells_in_basin[s];
730  for (int b = 0; b < m_n_basins; b++) {
731  if (n_shelf_cells_to_split[s * m_n_basins + b] > 0) {
732  n_shelf_numbers_to_add += 1;
733  add_shelf_instance[s * m_n_basins + b] = n_shelves + n_shelf_numbers_to_add;
734  m_log->message(3, "\nPICO, split ice shelf s=%d with bmax=%d "
735  "and b=%d and n=%d and si=%d\n", s, b0, b,
736  n_shelf_cells_to_split[s * m_n_basins + b],
737  add_shelf_instance[s * m_n_basins + b]);
738  }
739  }
740  }
741 
742  for (Points p(*m_grid); p; p.next()) {
743  const int i = p.i(), j = p.j();
744  if (cell_type.as_int(i, j) == MASK_FLOATING) {
745  int b = basin_mask.as_int(i, j);
746  int s = shelf_mask.as_int(i, j);
747  if (add_shelf_instance[s * m_n_basins + b] > 0) {
748  m_tmp(i, j) = add_shelf_instance[s * m_n_basins + b];
749  }
750  }
751  }
752 
753  shelf_mask.copy_from(m_tmp);
754 }
755 
756 /*!
757  * Compute distance to the grounding line.
758  */
760  const IceModelVec2Int &ice_rises,
761  bool exclude_ice_rises,
762  IceModelVec2Int &result) {
763 
764  IceModelVec::AccessList list{ &ice_rises, &ocean_mask, &result };
765 
766  result.set(-1);
767 
768  // Find the grounding line and the ice front and set result to 1 if ice shelf cell is
769  // next to the grounding line, Ice holes within the shelf are treated like ice shelf
770  // cells, if exclude_ice_rises is set then ice rises are also treated as ice shelf cells.
771 
772  ParallelSection loop(m_grid->com);
773  try {
774  for (Points p(*m_grid); p; p.next()) {
775  const int i = p.i(), j = p.j();
776 
777  if (ice_rises.as_int(i, j) == FLOATING or
778  ocean_mask.as_int(i, j) == 1 or
779  (exclude_ice_rises and ice_rises.as_int(i, j) == RISE)) {
780  // if this is an ice shelf cell (or an ice rise) or a hole in an ice shelf
781 
782  // label the shelf cells adjacent to the grounding line with 1
783  bool neighbor_to_land =
784  (ice_rises(i, j + 1) == CONTINENTAL or
785  ice_rises(i, j - 1) == CONTINENTAL or
786  ice_rises(i + 1, j) == CONTINENTAL or
787  ice_rises(i - 1, j) == CONTINENTAL or
788  ice_rises(i + 1, j + 1) == CONTINENTAL or
789  ice_rises(i + 1, j - 1) == CONTINENTAL or
790  ice_rises(i - 1, j + 1) == CONTINENTAL or
791  ice_rises(i - 1, j - 1) == CONTINENTAL);
792 
793  if (neighbor_to_land) {
794  // i.e. there is a grounded neighboring cell (which is not an ice rise!)
795  result(i, j) = 1;
796  } else {
797  result(i, j) = 0;
798  }
799  }
800  }
801  } catch (...) {
802  loop.failed();
803  }
804  loop.check();
805 
806  result.update_ghosts();
807 
808  eikonal_equation(result);
809 }
810 
811 /*!
812  * Compute distance to the calving front.
813  */
815  const IceModelVec2Int &ice_rises,
816  bool exclude_ice_rises,
817  IceModelVec2Int &result) {
818 
819  IceModelVec::AccessList list{ &ice_rises, &ocean_mask, &result };
820 
821  result.set(-1);
822 
823  ParallelSection loop(m_grid->com);
824  try {
825  for (Points p(*m_grid); p; p.next()) {
826  const int i = p.i(), j = p.j();
827 
828  if (ice_rises.as_int(i, j) == FLOATING or
829  ocean_mask.as_int(i, j) == 1 or
830  (exclude_ice_rises and ice_rises.as_int(i, j) == RISE)) {
831  // if this is an ice shelf cell (or an ice rise) or a hole in an ice shelf
832 
833  // label the shelf cells adjacent to the ice front with 1
834  auto M = ocean_mask.star(i, j);
835 
836  if (M.n == 2 or M.e == 2 or M.s == 2 or M.w == 2) {
837  // i.e. there is a neighboring open ocean cell
838  result(i, j) = 1;
839  } else {
840  result(i, j) = 0;
841  }
842  }
843  }
844  } catch (...) {
845  loop.failed();
846  }
847  loop.check();
848 
849  result.update_ghosts();
850 
851  eikonal_equation(result);
852 }
853 
854 /*!
855  * Find an approximate solution of the Eikonal equation on a given domain.
856  *
857  * To specify the problem, the input field (mask) should be filled with
858  *
859  * - negative values outside the domain
860  * - zeros within the domain
861  * - ones at "wave front" locations
862  *
863  * For example, to compute distances from the grounding line within ice shelves, fill
864  * generic ice shelf locations with zeros, set neighbors of the grounding line to 1, and
865  * the rest of the grid with -1 or some other negative number.
866  *
867  * Note that this implementation updates ghosts *every* iteration. We could speed this
868  * up by checking if a point at a boundary of the processor sub-domain was updated and
869  * update ghosts in those cases only.
870  *
871  * FIXME: replace this with a better algorithm.
872  */
874 
875  assert(mask.stencil_width() > 0);
876 
877  IceGrid::ConstPtr grid = mask.grid();
878 
879  double current_label = 1;
880  double continue_loop = 1;
881  while (continue_loop != 0) {
882 
883  continue_loop = 0;
884 
885  for (Points p(*grid); p; p.next()) {
886  const int i = p.i(), j = p.j();
887 
888  if (mask.as_int(i, j) == 0) {
889 
890  auto R = mask.star(i, j);
891 
892  if (R.ij == 0 and
893  (R.n == current_label or R.s == current_label or
894  R.e == current_label or R.w == current_label)) {
895  // i.e. this is an shelf cell with no distance assigned yet and with a neighbor
896  // that has a distance assigned
897  mask(i, j) = current_label + 1;
898  continue_loop = 1;
899  }
900  }
901  } // loop over grid points
902 
903  current_label++;
904  mask.update_ghosts();
905 
906  continue_loop = GlobalMax(grid->com, continue_loop);
907  }
908 }
909 
910 /*!
911  * Compute the mask identifying ice shelf "boxes" using distances to the grounding line
912  * and the calving front.
913  */
915  const IceModelVec2Int &shelf_mask, int max_number_of_boxes,
916  IceModelVec2Int &result) {
917 
918  IceModelVec::AccessList list{ &D_gl, &D_cf, &shelf_mask, &result };
919 
920  int n_shelves = static_cast<int>(shelf_mask.range()[1]) + 1;
921 
922  std::vector<double> GL_distance_max(n_shelves, 0.0);
923  std::vector<double> GL_distance_max1(n_shelves, 0.0);
924  std::vector<double> CF_distance_max(n_shelves, 0.0);
925  std::vector<double> CF_distance_max1(n_shelves, 0.0);
926 
927  for (Points p(*m_grid); p; p.next()) {
928  const int i = p.i(), j = p.j();
929 
930  int shelf_id = shelf_mask.as_int(i, j);
931  assert(shelf_id >= 0);
932  assert(shelf_id < n_shelves + 1);
933 
934  if (shelf_id == 0) {
935  // not at a shelf; skip to the next grid point
936  continue;
937  }
938 
939  if (D_gl(i, j) > GL_distance_max[shelf_id]) {
940  GL_distance_max[shelf_id] = D_gl(i, j);
941  }
942 
943  if (D_cf(i, j) > CF_distance_max[shelf_id]) {
944  CF_distance_max[shelf_id] = D_cf(i, j);
945  }
946  }
947 
948  // compute global maximums
949  GlobalMax(m_grid->com, GL_distance_max.data(), GL_distance_max1.data(), n_shelves);
950  GlobalMax(m_grid->com, CF_distance_max.data(), CF_distance_max1.data(), n_shelves);
951  // copy data
952  GL_distance_max = GL_distance_max1;
953  CF_distance_max = CF_distance_max1;
954 
955  double GL_distance_ref = *std::max_element(GL_distance_max.begin(), GL_distance_max.end());
956 
957  // compute the number of boxes in each shelf
958 
959  std::vector<int> n_boxes(n_shelves, 0);
960  const int n_min = 1;
961  const double zeta = 0.5;
962 
963  for (int k = 0; k < n_shelves; ++k) {
964  n_boxes[k] = n_min + round(pow((GL_distance_max[k] / GL_distance_ref), zeta) * (max_number_of_boxes - n_min));
965 
966  n_boxes[k] = std::min(n_boxes[k], max_number_of_boxes);
967  }
968 
969  result.set(0.0);
970 
971  for (Points p(*m_grid); p; p.next()) {
972  const int i = p.i(), j = p.j();
973 
974  int d_gl = D_gl.as_int(i, j);
975  int d_cf = D_cf.as_int(i, j);
976 
977  if (shelf_mask.as_int(i, j) > 0 and d_gl > 0.0 and d_cf > 0.0) {
978  int shelf_id = shelf_mask.as_int(i, j);
979  int n = n_boxes[shelf_id];
980 
981  // relative position on the shelf (ranges from 0 to 1), increasing towards the
982  // calving front (see equation 10 in the PICO paper)
983  double r = d_gl / (double)(d_gl + d_cf);
984 
985  double C = pow(1.0 - r, 2.0);
986  for (int k = 0; k < n; ++k) {
987  if ((n - k - 1) / (double)n <= C and C <= (n - k) / (double)n) {
988  result(i, j) = std::min(d_gl, k + 1);
989  }
990  }
991  }
992  }
993 }
994 
995 } // end of namespace ocean
996 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:101
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool next_to_ice(int i, int j) const
Ice-free margin (at least one of four neighbors has ice).
bool icy(int i, int j) const
bool ice_free_ocean(int i, int j) const
bool ocean(int i, int j) const
bool grounded(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
int as_int(int i, int j) const
stencils::Star< int > star(int i, int j) const
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
void copy_from(const IceModelVec2S &source)
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
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 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 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.
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void put_on_proc0(petsc::Vec &onp0) const
Puts a local IceModelVec2S on processor 0.
std::array< double, 2 > range() const
Result: min <- min(v[j]), max <- max(v[j]).
Definition: iceModelVec.cc:193
std::shared_ptr< petsc::Vec > allocate_proc0_copy() const
Definition: iceModelVec.cc:993
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
void compute_box_mask(const IceModelVec2Int &D_gl, const IceModelVec2Int &D_cf, const IceModelVec2Int &shelf_mask, int max_number_of_boxes, IceModelVec2Int &result)
std::shared_ptr< petsc::Vec > m_tmp_p0
void compute_ice_shelf_mask(const IceModelVec2Int &ice_rise_mask, const IceModelVec2Int &lake_mask, IceModelVec2Int &result)
const IceModelVec2Int & continental_shelf_mask() const
Definition: PicoGeometry.cc:58
void split_ice_shelves(const IceModelVec2CellType &cell_type, const IceModelVec2Int &basin_mask, const std::map< int, std::set< int > > &basin_neighbors, const std::vector< int > &most_shelf_cells_in_basin, const std::vector< int > &cfs_in_basins_per_shelf, int n_shelves, IceModelVec2Int &shelf_mask)
void compute_continental_shelf_mask(const IceModelVec2S &bed_elevation, const IceModelVec2Int &ice_rise_mask, double bed_elevation_threshold, IceModelVec2Int &result)
IceModelVec2Int m_distance_gl
const IceModelVec2Int & basin_mask() const
Definition: PicoGeometry.cc:74
IceModelVec2Int m_continental_shelf
void compute_ice_rises(const IceModelVec2CellType &cell_type, bool exclude_ice_rises, IceModelVec2Int &result)
IceModelVec2Int m_ice_rises
PicoGeometry(IceGrid::ConstPtr grid)
Definition: PicoGeometry.cc:32
void identify_calving_front_connection(const IceModelVec2CellType &cell_type, const IceModelVec2Int &basin_mask, const IceModelVec2Int &shelf_mask, int n_shelves, std::vector< int > &most_shelf_cells_in_basin, std::vector< int > &cfs_in_basins_per_shelf)
void compute_ocean_mask(const IceModelVec2CellType &cell_type, IceModelVec2Int &result)
void compute_distances_cf(const IceModelVec2Int &ocean_mask, const IceModelVec2Int &ice_rises, bool exclude_ice_rises, IceModelVec2Int &result)
const IceModelVec2Int & ice_shelf_mask() const
Definition: PicoGeometry.cc:66
void compute_distances_gl(const IceModelVec2Int &ocean_mask, const IceModelVec2Int &ice_rises, bool exclude_ice_rises, IceModelVec2Int &result)
void update(const IceModelVec2S &bed_elevation, const IceModelVec2CellType &cell_type)
Definition: PicoGeometry.cc:93
IceModelVec2Int m_lake_mask
IceModelVec2Int m_ocean_mask
std::map< int, std::set< int > > m_basin_neighbors
void compute_lakes(const IceModelVec2CellType &cell_type, IceModelVec2Int &result)
const IceModelVec2Int & box_mask() const
Definition: PicoGeometry.cc:62
IceModelVec2Int m_ice_shelves
IceModelVec2Int m_basin_mask
const IceModelVec2Int & ice_rise_mask() const
Definition: PicoGeometry.cc:70
std::map< int, std::set< int > > basin_neighbors(const IceModelVec2CellType &cell_type, const IceModelVec2Int &basin_mask)
IceModelVec2Int m_distance_cf
double * get()
Definition: Vec.cc:54
Wrapper around VecGetArray and VecRestoreArray.
Definition: Vec.hh:44
void label_connected_components(double *image, int nrows, int ncols, bool identify_icebergs, int mask_grounded, int first_label)
#define PISM_ERROR_LOCATION
static const double b0
Definition: exactTestL.cc:41
#define n
Definition: exactTestM.c:37
bool ice_free_ocean(int M)
Definition: Mask.hh:60
bool ocean(int M)
An ocean cell (floating ice or ice-free).
Definition: Mask.hh:39
static void relabel(RelabelingType type, double threshold, IceModelVec2Int &mask)
void eikonal_equation(IceModelVec2Int &mask)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
std::string printf(const char *format,...)
static const double k
Definition: exactTestP.cc:45
@ MASK_FLOATING
Definition: Mask.hh:33
@ MASK_ICE_FREE_OCEAN
Definition: Mask.hh:34
@ CRITICAL
Definition: IO_Flags.hh:70
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
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)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
bool grid_edge(const IceGrid &grid, int i, int j)
Definition: IceGrid.hh:350
const int C[]
Definition: ssafd_code.cc:44
std::string filename
Definition: options.hh:33