PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
IcebergRemoverFEM.cc
Go to the documentation of this file.
1 /* Copyright (C) 2021, 2022, 2023 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 <cassert>
21 
22 #include "pism/util/petscwrappers/DM.hh"
23 #include "pism/util/petscwrappers/Vec.hh"
24 #include "pism/util/connected_components.hh"
25 
26 #include "pism/frontretreat/util/IcebergRemoverFEM.hh"
27 
28 #include "pism/util/fem/Element.hh"
29 #include "pism/util/array/CellType.hh"
30 #include "pism/util/Mask.hh"
31 #include "pism/util/interpolation.hh"
32 
33 namespace pism {
34 namespace calving {
35 
36 IcebergRemoverFEM::IcebergRemoverFEM(std::shared_ptr<const Grid> grid)
37  : IcebergRemover(grid),
38  m_mask(grid, "temporary_mask") {
40 }
41 
42 /*! Remove "icebergs" using the finite element notion of connectivity: two elements are
43  * connected if they share a boundary.
44  *
45  * 1. Loop over elements and create a mask that will be used to determine connectivity
46  * between elements.
47  *
48  * - an element is a "grounded ice element" if all nodes are icy and are either grounded
49  * or belong to the set of Dirichlet nodes
50  *
51  * - an element is "floating ice" if all nodes are icy and at least one node is
52  * "floating ice"
53  *
54  * - all other elements are ice-free
55  *
56  * 2. Label connected components, identifying "icebergs".
57  *
58  * Once "iceberg" elements are labeled we need to remove *nodes* that belong to icebergs
59  * but *do not* belong to any elements connected to grounded ice.
60  *
61  * 3. Create a mask filled with zeros. Loop over elements and add 1 to nodes of all
62  * "iceberg" elements. Add -1 to all nodes of "grounded" elements.
63  *
64  * 4. Now loop over all nodes and remove nodes with positive mask values.
65  *
66  */
68  array::CellType1 &cell_type,
69  array::Scalar &ice_thickness) {
70  const int
71  mask_grounded_ice = 1,
72  mask_floating_ice = 2;
73 
74  int bc_mask_nodal[fem::q1::n_chi];
75  int cell_type_nodal[fem::q1::n_chi];
76 
77  assert(bc_mask.stencil_width() >= 1);
78  assert(cell_type.stencil_width() >= 1);
79 
80  array::AccessScope list{&bc_mask, &cell_type, &m_iceberg_mask};
81 
83 
84  // prepare the iceberg mask: the value at (i, j) describes an *element* with (i,j) as a
85  // lower left corner
86  {
87  // loop over all nodes in a local sub-domain
88  for (auto p = m_grid->points(); p; p.next()) {
89  const int i = p.i(), j = p.j();
90 
91  element.reset(i, j);
92  // the following two calls use ghost values
93  element.nodal_values(bc_mask, bc_mask_nodal);
94  element.nodal_values(cell_type, cell_type_nodal);
95 
96  // check if all nodes are icy
97  bool icy = true;
98  for (int n = 0; icy and n < fem::q1::n_chi; ++n) {
99  icy &= mask::icy(cell_type_nodal[n]);
100  }
101 
102  if (icy) {
103  // This is an icy element: check if all nodes are grounded or are a part of the
104  // set of Dirichlet nodes
105  bool grounded = true;
106  for (int n = 0; grounded and n < fem::q1::n_chi; ++n) {
107  grounded &= (mask::grounded(cell_type_nodal[n]) or bc_mask_nodal[n] == 1);
108  }
109 
110  m_iceberg_mask(i, j) = grounded ? mask_grounded_ice : mask_floating_ice;
111  } else {
112  // This is an ice-free element.
113  m_iceberg_mask(i, j) = 0;
114  }
115  } // end of the loop over local nodes
116  } // end of the block preparing the mask
117 
118  // Identify icebergs using serial code on processor 0:
119  {
121 
122  ParallelSection rank0(m_grid->com);
123  try {
124  if (m_grid->rank() == 0) {
125  petsc::VecArray mask_p0(*m_mask_p0);
126  label_connected_components(mask_p0.get(), m_grid->My(), m_grid->Mx(),
127  true, mask_grounded_ice);
128  }
129  } catch (...) {
130  rank0.failed();
131  }
132  rank0.check();
133 
135  // note: this will update ghosts of m_iceberg_mask
136  }
137 
138  // create a mask indicating if a *node* should be removed
139  {
140  DMDALocalInfo info;
141  {
142  auto da = m_grid->get_dm(1, 0); // dof = 1, stencil_width = 0
143  PetscErrorCode ierr = DMDAGetLocalInfo(*da, &info);
144  if (ierr != 0) {
145  throw std::runtime_error("Failed to get DMDA info");
146  }
147  }
148 
149  m_mask.set(0);
150  list.add(m_mask);
151  double **M = m_mask.array();
152 
153  double mask_iceberg[] = {1.0, 1.0, 1.0, 1.0};
154  double mask_grounded[] = {-1.0, -1.0, -1.0, -1.0};
155 
156  // loop over all the elements that have at least one owned node
157  for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
158  for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
159  element.reset(i, j);
160 
161  // the following two calls use ghost values
162  element.nodal_values(bc_mask, bc_mask_nodal);
163  element.nodal_values(cell_type, cell_type_nodal);
164 
165  // check if all nodes are icy
166  bool icy = true;
167  for (int n = 0; icy and n < fem::q1::n_chi; ++n) {
168  icy &= mask::icy(cell_type_nodal[n]);
169  }
170 
171  if (icy) {
172  // check if all nodes are grounded or are a part of the set of Dirichlet nodes
173  bool grounded = true;
174  for (int n = 0; grounded and n < fem::q1::n_chi; ++n) {
175  grounded &= (mask::grounded(cell_type_nodal[n]) or bc_mask_nodal[n] == 1);
176  }
177 
178  if (m_iceberg_mask(i, j) == 1) {
179  // this is an iceberg element
180  element.add_contribution(mask_iceberg, M);
181  } else {
182  element.add_contribution(mask_grounded, M);
183  }
184  }
185  }
186  } // end of the loop over elements
187  } // end of the block identifying nodes to remove
188 
189  // loop over all *nodes* and modify ice thickness and mask
190  {
191  list.add(ice_thickness);
192 
193  for (auto p = m_grid->points(); p; p.next()) {
194  const int i = p.i(), j = p.j();
195 
196  if (m_mask(i, j) > 0) {
197  ice_thickness(i,j) = 0.0;
198  cell_type(i,j) = MASK_ICE_FREE_OCEAN;
199  }
200  }
201  }
202 
203  // update ghosts of the mask and the ice thickness (then surface
204  // elevation can be updated redundantly)
205  cell_type.update_ghosts();
206  ice_thickness.update_ghosts();
207 }
208 
209 } // end of namespace calving
210 } // end of namespace pism
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
void failed()
Indicates a failure of a parallel section.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void set_interpolation_type(InterpolationType type)
Definition: Array.cc:179
void get_from_proc0(petsc::Vec &onp0)
Gets a local Array2 from processor 0.
Definition: Array.cc:1095
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void put_on_proc0(petsc::Vec &onp0) const
Puts a local array::Scalar on processor 0.
Definition: Array.cc:1052
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
unsigned int stencil_width() const
Get the stencil width of the current Array. Returns 0 if ghosts are not available.
Definition: Array.cc:331
IcebergRemoverFEM(std::shared_ptr< const Grid > g)
void update_impl(const array::Scalar &bc_mask, array::CellType1 &cell_type, array::Scalar &ice_thickness)
std::shared_ptr< petsc::Vec > m_mask_p0
PISM iceberg remover.
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
Definition: Element.cc:187
void add_contribution(const T *local, T **y_global) const
Add the values of element-local contributions y to the global vector y_global.
Definition: Element.hh:231
void nodal_values(const array::Scalar &x_global, int *result) const
Get nodal values of an integer mask.
Definition: Element.cc:176
Q1 element with sides parallel to X and Y axes.
Definition: Element.hh:257
The 1-point Gaussian quadrature on the square [-1,1]*[-1,1].
Definition: Quadrature.hh:89
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 n
Definition: exactTestM.c:37
const int n_chi
Definition: FEM.hh:191
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:48
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:44
@ MASK_ICE_FREE_OCEAN
Definition: Mask.hh:35