PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
node_types.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 2018, 2020, 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 "pism/util/node_types.hh"
21 
22 #include "pism/util/array/Scalar.hh"
23 #include "pism/util/array/Scalar.hh"
24 #include "pism/util/Grid.hh"
25 #include "pism/util/error_handling.hh"
26 
27 namespace pism {
28 
29 /**
30  Identify node types of all the nodes (grid points).
31 
32  Uses width=1 ghosts of `ice_thickness` (a "box" stencil).
33 
34  A node is can be *interior*, *boundary*, or *exterior*.
35 
36  An element is considered *icy* if at least three of its nodes have ice thickness above
37  `thickness_threshold`.
38 
39  A node is considered *interior* if all of the elements it belongs to are icy.
40 
41  A node is considered *boundary* if it is not interior and at least one element it belongs to is
42  icy.
43 
44  A node is considered *exterior* if it is neither interior nor boundary.
45 
46  Now an element "face" (side) is a part of a boundary if and only if both of its nodes are
47  boundary nodes.
48 
49 Cell layout:
50 ~~~
51 (i-1,j+1) +-------N--------+ (i+1,j+1)
52  | | |
53  | NW | NE |
54  | | |
55  W-----(i,j)------E
56  | | |
57  | SW | SE |
58  | | |
59 (i-1,j-1) +-------S--------+ (i+1,j-1)
60 ~~~
61  */
62 void compute_node_types(const array::Scalar1 &ice_thickness,
63  double thickness_threshold,
64  array::Scalar &result) {
65 
66  auto grid = ice_thickness.grid();
67 
68  const double &H_min = thickness_threshold;
69 
70  array::AccessScope list{&ice_thickness, &result};
71 
72  ParallelSection loop(grid->com);
73  try {
74  for (auto p = grid->points(); p; p.next()) {
75  const int i = p.i(), j = p.j();
76 
77  auto H = ice_thickness.box(i, j);
78 
79  // flags indicating whether the current node and its neighbors are "icy"
81 
82  icy.c = static_cast<int>(H.c >= H_min);
83  icy.nw = static_cast<int>(H.nw >= H_min);
84  icy.n = static_cast<int>(H.n >= H_min);
85  icy.ne = static_cast<int>(H.ne >= H_min);
86  icy.e = static_cast<int>(H.e >= H_min);
87  icy.se = static_cast<int>(H.se >= H_min);
88  icy.s = static_cast<int>(H.s >= H_min);
89  icy.sw = static_cast<int>(H.sw >= H_min);
90  icy.w = static_cast<int>(H.w >= H_min);
91 
92  // flags indicating whether neighboring elements are "icy" (an element is icy if at
93  // least three of its nodes are icy)
94  const bool
95  ne_element_is_icy = (icy.c + icy.e + icy.ne + icy.n) >= 3,
96  nw_element_is_icy = (icy.c + icy.n + icy.nw + icy.w) >= 3,
97  sw_element_is_icy = (icy.c + icy.w + icy.sw + icy.s) >= 3,
98  se_element_is_icy = (icy.c + icy.s + icy.se + icy.e) >= 3;
99 
100  if (ne_element_is_icy and nw_element_is_icy and
101  sw_element_is_icy and se_element_is_icy) {
102  // all four elements are icy: we are at an interior node
103  result(i, j) = NODE_INTERIOR;
104  } else if (icy.c != 0) {
105  // the current node is icy: we are at a boundary
106  result(i, j) = NODE_BOUNDARY;
107  } else {
108  // all elements are ice-free: we are at an exterior node
109  result(i, j) = NODE_EXTERIOR;
110  }
111 
112  }
113  } catch (...) {
114  loop.failed();
115  }
116  loop.check();
117 
118  result.update_ghosts();
119 }
120 
121 } // end of namespace pism
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
stencils::Box< T > box(int i, int j) const
Definition: Array2D.hh:93
std::shared_ptr< const Grid > grid() const
Definition: Array.cc:132
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:48
@ NODE_BOUNDARY
Definition: node_types.hh:35
@ NODE_EXTERIOR
Definition: node_types.hh:36
@ NODE_INTERIOR
Definition: node_types.hh:34
void compute_node_types(const array::Scalar1 &ice_thickness, double thickness_threshold, array::Scalar &result)
Definition: node_types.cc:62