PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
IcebergRemoverFEM.cc
Go to the documentation of this file.
1 /* Copyright (C) 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 <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 "IcebergRemoverFEM.hh"
27 
28 #include "pism/util/fem/Element.hh"
29 #include "pism/util/IceModelVec2CellType.hh"
30 #include "pism/util/Mask.hh"
31 
32 namespace pism {
33 namespace calving {
34 
36  : IcebergRemover(grid),
37  m_mask(grid, "temporary_mask", WITHOUT_GHOSTS) {
38  // empty
39 }
40 
41 /*! Remove "icebergs" using the finite element notion of connectivity: two elements are
42  * connected if they share a boundary.
43  *
44  * 1. Loop over elements and create a mask that will be used to determine connectivity
45  * between elements.
46  *
47  * - an element is a "grounded ice element" if all nodes are icy and are either grounded
48  * or belong to the set of Dirichlet nodes
49  *
50  * - an element is "floating ice" if all nodes are icy and at least one node is
51  * "floating ice"
52  *
53  * - all other elements are ice-free
54  *
55  * 2. Label connected components, identifying "icebergs".
56  *
57  * Once "iceberg" elements are labeled we need to remove *nodes* that belong to icebergs
58  * but *do not* belong to any elements connected to grounded ice.
59  *
60  * 3. Create a mask filled with zeros. Loop over elements and add 1 to nodes of all
61  * "iceberg" elements. Add -1 to all nodes of "grounded" elements.
62  *
63  * 4. Now loop over all nodes and remove nodes with positive mask values.
64  *
65  */
67  IceModelVec2CellType &pism_mask,
68  IceModelVec2S &ice_thickness) {
69  const int
70  mask_grounded_ice = 1,
71  mask_floating_ice = 2;
72 
73  int bc_mask_nodal[fem::q1::n_chi];
74  int pism_mask_nodal[fem::q1::n_chi];
75 
76  assert(bc_mask.stencil_width() >= 1);
77  assert(pism_mask.stencil_width() >= 1);
78 
79  IceModelVec::AccessList list{&bc_mask, &pism_mask, &m_iceberg_mask};
80 
82 
83  // prepare the iceberg mask: the value at (i, j) describes an *element* with (i,j) as a
84  // lower left corner
85  {
86  // loop over all nodes in a local sub-domain
87  for (Points p(*m_grid); p; p.next()) {
88  const int i = p.i(), j = p.j();
89 
90  element.reset(i, j);
91  // the following two calls use ghost values
92  element.nodal_values(bc_mask, bc_mask_nodal);
93  element.nodal_values(pism_mask, pism_mask_nodal);
94 
95  // check if all nodes are icy
96  bool icy = true;
97  for (int n = 0; icy and n < fem::q1::n_chi; ++n) {
98  icy &= mask::icy(pism_mask_nodal[n]);
99  }
100 
101  if (icy) {
102  // This is an icy element: check if all nodes are grounded or are a part of the
103  // set of Dirichlet nodes
104  bool grounded = true;
105  for (int n = 0; grounded and n < fem::q1::n_chi; ++n) {
106  grounded &= (mask::grounded(pism_mask_nodal[n]) or bc_mask_nodal[n] == 1);
107  }
108 
109  m_iceberg_mask(i, j) = grounded ? mask_grounded_ice : mask_floating_ice;
110  } else {
111  // This is an ice-free element.
112  m_iceberg_mask(i, j) = 0;
113  }
114  } // end of the loop over local nodes
115  } // end of the block preparing the mask
116 
117  // Identify icebergs using serial code on processor 0:
118  {
120 
121  ParallelSection rank0(m_grid->com);
122  try {
123  if (m_grid->rank() == 0) {
124  petsc::VecArray mask_p0(*m_mask_p0);
125  label_connected_components(mask_p0.get(), m_grid->My(), m_grid->Mx(),
126  true, mask_grounded_ice);
127  }
128  } catch (...) {
129  rank0.failed();
130  }
131  rank0.check();
132 
134  // note: this will update ghosts of m_iceberg_mask
135  }
136 
137  // create a mask indicating if a *node* should be removed
138  {
139  DMDALocalInfo info;
140  {
141  auto da = m_grid->get_dm(1, 0); // dof = 1, stencil_width = 0
142  PetscErrorCode ierr = DMDAGetLocalInfo(*da, &info);
143  if (ierr != 0) {
144  throw std::runtime_error("Failed to get DMDA info");
145  }
146  }
147 
148  m_mask.set(0);
149  list.add(m_mask);
150  double **M = m_mask.array();
151 
152  double mask_iceberg[] = {1.0, 1.0, 1.0, 1.0};
153  double mask_grounded[] = {-1.0, -1.0, -1.0, -1.0};
154 
155  // loop over all the elements that have at least one owned node
156  for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
157  for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
158  element.reset(i, j);
159 
160  // the following two calls use ghost values
161  element.nodal_values(bc_mask, bc_mask_nodal);
162  element.nodal_values(pism_mask, pism_mask_nodal);
163 
164  // check if all nodes are icy
165  bool icy = true;
166  for (int n = 0; icy and n < fem::q1::n_chi; ++n) {
167  icy &= mask::icy(pism_mask_nodal[n]);
168  }
169 
170  if (icy) {
171  // check if all nodes are grounded or are a part of the set of Dirichlet nodes
172  bool grounded = true;
173  for (int n = 0; grounded and n < fem::q1::n_chi; ++n) {
174  grounded &= (mask::grounded(pism_mask_nodal[n]) or bc_mask_nodal[n] == 1);
175  }
176 
177  if (m_iceberg_mask(i, j) == 1) {
178  // this is an iceberg element
179  element.add_contribution(mask_iceberg, M);
180  } else {
181  element.add_contribution(mask_grounded, M);
182  }
183  }
184  }
185  } // end of the loop over elements
186  } // end of the block identifying nodes to remove
187 
188  // loop over all *nodes* and modify ice thickness and mask
189  {
190  list.add(ice_thickness);
191 
192  for (Points p(*m_grid); p; p.next()) {
193  const int i = p.i(), j = p.j();
194 
195  if (m_mask(i, j) > 0) {
196  ice_thickness(i,j) = 0.0;
197  pism_mask(i,j) = MASK_ICE_FREE_OCEAN;
198  }
199  }
200  }
201 
202  // update ghosts of the mask and the ice thickness (then surface
203  // elevation can be updated redundantly)
204  pism_mask.update_ghosts();
205  ice_thickness.update_ghosts();
206 }
207 
208 } // end of namespace calving
209 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
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
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.
void put_on_proc0(petsc::Vec &onp0) const
Puts a local IceModelVec2S on processor 0.
void failed()
Indicates a failure of a parallel section.
IcebergRemoverFEM(IceGrid::ConstPtr g)
void update_impl(const IceModelVec2Int &bc_mask, IceModelVec2CellType &pism_mask, IceModelVec2S &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 IceModelVec2Int &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:47
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:43
@ MASK_ICE_FREE_OCEAN
Definition: Mask.hh:34
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49