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