tIcebergRemover.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
---
tIcebergRemover.cc (3670B)
---
1 /* Copyright (C) 2013, 2014, 2015, 2016, 2017, 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 "IcebergRemover.hh"
21 #include "pism/util/connected_components.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/IceGrid.hh"
26 #include "pism/util/IceModelVec2CellType.hh"
27
28 namespace pism {
29 namespace calving {
30
31 IcebergRemover::IcebergRemover(IceGrid::ConstPtr g)
32 : Component(g),
33 m_iceberg_mask(m_grid, "iceberg_mask", WITHOUT_GHOSTS){
34
35 m_mask_p0 = m_iceberg_mask.allocate_proc0_copy();
36 }
37
38 IcebergRemover::~IcebergRemover() {
39 // empty
40 }
41
42 void IcebergRemover::init() {
43 }
44
45 /**
46 * Use PISM's ice cover mask to update ice thickness, removing "icebergs".
47 *
48 * @param[in,out] pism_mask PISM's ice cover mask
49 * @param[in,out] ice_thickness ice thickness
50 */
51 void IcebergRemover::update(const IceModelVec2Int &bc_mask,
52 IceModelVec2CellType &mask,
53 IceModelVec2S &ice_thickness) {
54 const int
55 mask_grounded_ice = 1,
56 mask_floating_ice = 2;
57
58 // prepare the mask that will be handed to the connected component
59 // labeling code:
60 {
61 m_iceberg_mask.set(0.0);
62
63 IceModelVec::AccessList list{&mask, &m_iceberg_mask, &bc_mask};
64
65 for (Points p(*m_grid); p; p.next()) {
66 const int i = p.i(), j = p.j();
67
68 if (mask.grounded_ice(i,j) == true) {
69 m_iceberg_mask(i,j) = mask_grounded_ice;
70 } else if (mask.floating_ice(i,j) == true) {
71 m_iceberg_mask(i,j) = mask_floating_ice;
72 }
73 }
74
75 // Mark icy Dirichlet B.C. cells as "grounded" because we don't want them removed.
76 for (Points p(*m_grid); p; p.next()) {
77 const int i = p.i(), j = p.j();
78
79 if (bc_mask(i, j) > 0.5 and mask.icy(i, j)) {
80 m_iceberg_mask(i, j) = mask_grounded_ice;
81 }
82 }
83 }
84
85 // identify icebergs using serial code on processor 0:
86 {
87 m_iceberg_mask.put_on_proc0(*m_mask_p0);
88
89 ParallelSection rank0(m_grid->com);
90 try {
91 if (m_grid->rank() == 0) {
92 petsc::VecArray mask_p0(*m_mask_p0);
93 label_connected_components(mask_p0.get(), m_grid->My(), m_grid->Mx(), true, mask_grounded_ice);
94 }
95 } catch (...) {
96 rank0.failed();
97 }
98 rank0.check();
99
100 m_iceberg_mask.get_from_proc0(*m_mask_p0);
101 }
102
103 // correct ice thickness and the cell type mask using the resulting
104 // "iceberg" mask:
105 {
106 IceModelVec::AccessList list{&ice_thickness, &mask, &m_iceberg_mask, &bc_mask};
107
108 for (Points p(*m_grid); p; p.next()) {
109 const int i = p.i(), j = p.j();
110
111 if (m_iceberg_mask(i,j) > 0.5 && bc_mask(i,j) < 0.5) {
112 ice_thickness(i,j) = 0.0;
113 mask(i,j) = MASK_ICE_FREE_OCEAN;
114 }
115 }
116 }
117
118 // update ghosts of the mask and the ice thickness (then surface
119 // elevation can be updated redundantly)
120 mask.update_ghosts();
121 ice_thickness.update_ghosts();
122 }
123
124 } // end of namespace calving
125 } // end of namespace pism