tPicoGeometry.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tPicoGeometry.cc (19955B)
---
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
30 PicoGeometry::PicoGeometry(IceGrid::ConstPtr grid)
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",
45 {OCEAN, RISE, CONTINENTAL, FLOATING});
46 m_ice_rises.metadata().set_string("flag_meanings",
47 "ocean ice_rise continental_ice_sheet, floating_ice");
48
49 m_tmp_p0 = m_tmp.allocate_proc0_copy();
50 }
51
52 PicoGeometry::~PicoGeometry() {
53 // empty
54 }
55
56 const IceModelVec2Int &PicoGeometry::continental_shelf_mask() const {
57 return m_continental_shelf;
58 }
59
60 const IceModelVec2Int &PicoGeometry::box_mask() const {
61 return m_boxes;
62 }
63
64 const IceModelVec2Int &PicoGeometry::ice_shelf_mask() const {
65 return m_ice_shelves;
66 }
67
68 const IceModelVec2Int &PicoGeometry::ice_rise_mask() const {
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 {
96 m_ice_rises.update_ghosts();
97 m_ocean_mask.update_ghosts();
98
99 compute_distances_gl(m_ocean_mask, m_ice_rises, exclude_ice_rises, m_distance_gl);
100
101 compute_distances_cf(m_ocean_mask, m_ice_rises, exclude_ice_rises, m_distance_cf);
102 }
103
104 // these two could be done at the same time
105 {
106 compute_ice_shelf_mask(m_ice_rises, m_lake_mask, m_ice_shelves);
107
108 compute_continental_shelf_mask(bed_elevation, m_ice_rises, continental_shelf_depth, m_continental_shelf);
109 }
110
111 compute_box_mask(m_distance_gl, m_distance_cf, m_ice_shelves, n_boxes, m_boxes);
112 }
113
114
115 enum RelabelingType {BY_AREA, AREA_THRESHOLD};
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 */
212 void PicoGeometry::label_tmp() {
213 m_tmp.put_on_proc0(*m_tmp_p0);
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
226 m_tmp.get_from_proc0(*m_tmp_p0);
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 */
243 void PicoGeometry::compute_lakes(const IceModelVec2CellType &cell_type, IceModelVec2Int &result) {
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 {
268 m_tmp.put_on_proc0(*m_tmp_p0);
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
281 m_tmp.get_from_proc0(*m_tmp_p0);
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
316 relabel(AREA_THRESHOLD,
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 */
342 void PicoGeometry::compute_continental_shelf_mask(const IceModelVec2S &bed_elevation,
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 {
364 m_tmp.put_on_proc0(*m_tmp_p0);
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
377 m_tmp.get_from_proc0(*m_tmp_p0);
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 */
449 void PicoGeometry::compute_ocean_mask(const IceModelVec2CellType &cell_type, IceModelVec2Int &result) {
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
470 void PicoGeometry::compute_distances_gl(const IceModelVec2Int &ocean_mask,
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
522 void PicoGeometry::compute_distances_cf(const IceModelVec2Int &ocean_mask,
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 */
579 void eikonal_equation(IceModelVec2Int &mask) {
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
616 void PicoGeometry::compute_box_mask(const IceModelVec2Int &D_gl, const IceModelVec2Int &D_cf,
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