PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
grounded_cell_fraction.cc
Go to the documentation of this file.
1 /* Copyright (C) 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
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 #include <cmath> // fabs, isnan
22
23 #include "pism/geometry/grounded_cell_fraction.hh"
24
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/pism_utilities.hh" // clip
27 #include "pism/util/array/Scalar.hh"
28
29 namespace pism {
30
31 struct Point {
32  double x;
33  double y;
34 };
35
36 /*!
37  * Consider the right triangle A(0,0) - B(1,0) - C(0,1).
38  *
39  * Define a linear function z = a + (b - a) * x + (c - a) * y, where a, b, and c are its
40  * values at points A, B, and C, respectively.
41  *
42  * Our goal is to find the fraction of the area of ABC in where z > 0.
43  *
44  * We note that z(x,y) is continuous, so unless a, b, and c have the same sign the line
45  *
46  * z = 0
47  *
48  * will intersect exactly two of the sides (possibly at a node of ABC).
49  *
50  * So, if the line (z = 0) does not intersect BC, for example, then it has to intersect AB
51  * and AC.
52  *
53  * This method can be applied to arbitrary triangles. It does not even matter if values at
54  * triangle nodes (a, b, c) are listed in the same order. (For any two triangles on a
55  * plane there exists an affine map that takes one to the other. Also, affine maps
56  * preserve ratios of areas of figures.)
57  */
58
59 /*!
60  * Compute the area of a triangle on a plane using the "shoelace formula."
61  */
62 static inline double triangle_area(const Point &a, const Point &b, const Point &c) {
63  // note: fabs should not be needed since we traverse all triangle nodes
64  // counter-clockwise, but it is good to be safe
65  return 0.5 * std::fabs((a.x - c.x) * (b.y - a.y) - (a.x - b.x) * (c.y - a.y));
66 }
67
68 /*!
69  * Compute the coordinates of the intersection of (z = 0) with the side AB.
70  */
71 Point intersect_ab(double a, double b) {
72  if (a != b) {
73  return {a / (a - b), 0.0};
74  }
75  return {-1.0, -1.0}; // no intersection
76 }
77
78 /*!
79  * Compute the coordinates of the intersection of (z = 0) with the side BC.
80  */
81 Point intersect_bc(double b, double c) {
82  if (b != c) {
83  return {c / (c - b), b / (b - c)};
84  }
85  return {-1.0, -1.0}; // no intersection
86 }
87
88 /*!
89  * Compute the coordinates of the intersection of (z = 0) with the side AC.
90  */
91 Point intersect_ac(double a, double c) {
92  if (a != c) {
93  return {0.0, a / (a - c)};
94  }
95  return {-1.0, -1.0}; // no intersection
96 }
97
98 /*!
99  * Return true if a point p is not a valid point on a side of the triangle
100  * (0,0)-(1,0)-(0,1).
101  *
102  * This is not a complete test (for example, it does not check if y = 1 - x for points on
103  * the hypotenuse). The point of this is to exclude the kinds of invalid points we are
104  * likely to see, not all of them.
105  *
106  * Note that we use (-1, -1) to indicate "invalid" points in the rest of the code and
107  * these are easy to detect: they require only one comparison.
108  */
109 bool invalid(const Point &p) {
110  return (p.x < 0.0 or p.x > 1.0 or p.y < 0.0 or p.y > 1.0);
111 }
112
113 /*!
114  * Return true if two points are the same.
115  */
116 static bool same(const Point &a, const Point &b) {
117  const double threshold = 1e-12;
118  return std::fabs(a.x - b.x) < threshold and std::fabs(a.y - b.y) < threshold;
119 }
120
121 /*!
122  * Consider the right triangle A(0,0) - B(1,0) - C(0,1).
123  *
124  * Define a linear function z = a + (b - a) * x + (c - a) * y, where a, b, and c are its
125  * values at points A, B, and C, respectively.
126  *
127  * Our goal is to find the fraction of the triangle ABC in where z > 0.
128  *
129  * This corresponds to the grounded area fraction if z is the flotation criterion
130  * function.
131  */
132 double grounded_area_fraction(double a, double b, double c) {
133
134  assert(std::isfinite(a));
135  assert(std::isfinite(b));
136  assert(std::isfinite(c));
137
138  if (a > 0.0 and b > 0.0 and c > 0.0) {
139  return 1.0;
140  }
141
142  if (a <= 0.0 and b <= 0.0 and c <= 0.0) {
143  return 0.0;
144  }
145
146  // the area of the triangle (0,0)-(1,0)-(0,1)
147  const double total_area = 0.5;
148
149  const Point
150  ab = intersect_ab(a, b),
151  bc = intersect_bc(b, c),
152  ac = intersect_ac(a, c);
153
154  if (invalid(bc)) {
155  assert(not (invalid(ab) or invalid(ac)));
156
157  double ratio = triangle_area({0.0, 0.0}, ab, ac) / total_area;
158  assert((ratio >= 0.0) and (ratio <= 1.0));
159
160  if (a > 0.0) {
161  return ratio;
162  }
163  return 1.0 - ratio;
164  }
165
166  if (invalid(ac)) {
167  assert(not (invalid(ab) or invalid(bc)));
168
169  double ratio = triangle_area({1.0, 0.0}, bc, ab) / total_area;
170  assert((ratio >= 0.0) and (ratio <= 1.0));
171
172  if (b > 0.0) {
173  return ratio;
174  }
175  return 1.0 - ratio;
176  }
177
178  if (invalid(ab)) {
179  assert(not (invalid(bc) or invalid(ac)));
180
181  double ratio = triangle_area({0.0, 1.0}, ac, bc) / total_area;
182  assert((ratio >= 0.0) and (ratio <= 1.0));
183
184  if (c > 0.0) {
185  return ratio;
186  }
187  return 1.0 - ratio;
188  }
189
190  // Note that we know that ab, bc, and ac are all valid.
191
192  // the a == 0 case, the line F = 0 goes through A
193  if (same(ab, ac)) {
194  double ratio = triangle_area({1.0, 0.0}, bc, ab) / total_area;
195  assert((ratio >= 0.0) and (ratio <= 1.0));
196
197  if (b > 0.0) {
198  return ratio;
199  }
200  return 1.0 - ratio;
201  }
202
203  // the b == 0 case and the c == 0 case
204  if (same(ab, bc) or same(ac, bc)) {
205  double ratio = triangle_area({0.0, 0.0}, ab, ac) / total_area;
206  assert((ratio >= 0.0) and (ratio <= 1.0));
207
208  if (a > 0.0) {
209  return ratio;
210  }
211  return 1.0 - ratio;
212  }
213
214  // Note: the case of F=0 coinciding with a side of the triangle is covered by if clauses
215  // above. For example, when F=0 coincides with AC, we have a = c = 0 and intersect_ac(a, c)
216  // returns an invalid intersection point.
217
219  "the logic in grounded_area_fraction failed! Please submit a bug report.");
220 }
221
222 /*!
223  * The flotation criterion.
224  */
225 static double F(double SL, double B, double H, double alpha) {
226  double
227  water_depth = SL - B,
228  shelf_depth = H * alpha;
229  return shelf_depth - water_depth;
230 }
231
232
233 /*!
234  * Compute the flotation criterion at all the points in the box stencil.
235  */
237 static Box F(const Box &SL, const Box &B, const Box &H, double alpha) {
238  return {F(SL.c, B.c, H.c, alpha),
239  F(SL.n, B.n, H.n, alpha),
240  F(SL.nw, B.nw, H.nw, alpha),
241  F(SL.w, B.w, H.w, alpha),
242  F(SL.sw, B.sw, H.sw, alpha),
243  F(SL.s, B.s, H.s, alpha),
244  F(SL.se, B.se, H.se, alpha),
245  F(SL.e, B.e, H.e, alpha),
246  F(SL.ne, B.ne, H.ne, alpha)};
247 }
248
249 /*!
250  * @param[in] ice_density ice density, kg/m3
251  * @param[in] ocean_density ocean_density, kg/m3
252  * @param[in] sea_level sea level (flotation) elevation, m
253  * @param[in] ice_thickness ice thickness, m
254  * @param[in] bed_topography bed elevation, m
255  * @param[out] result grounded cell fraction, between 0 (floating) and 1 (grounded)
256  */
257 void compute_grounded_cell_fraction(double ice_density,
258  double ocean_density,
259  const array::Scalar1 &sea_level,
260  const array::Scalar1 &ice_thickness,
261  const array::Scalar1 &bed_topography,
262  array::Scalar &result) {
263  auto grid = result.grid();
264  double alpha = ice_density / ocean_density;
265
266  array::AccessScope list{&sea_level, &ice_thickness, &bed_topography, &result};
267
268  ParallelSection loop(grid->com);
269  try {
270  for (auto p = grid->points(); p; p.next()) {
271  const int i = p.i(), j = p.j();
272
273  /*
274  NW----------------N----------------NE
275  | | |
276  | | |
277  | nw--------n--------ne |
278  | | | | |
279  | | | | |
280  W--------w--------o--------e--------E
281  | | | | |
282  | | | | |
283  | sw--------s--------se |
284  | | |
285  | | |
286  SW----------------S----------------SE
287  */
288
289  // compute the floatation function at 8 points surrounding the current grid point
291  {
292  auto S = sea_level.box(i, j);
293  auto H = ice_thickness.box(i, j);
294  auto B = bed_topography.box(i, j);
295
296  auto x = F(S, B, H, alpha);
297
298  f.c = x.c;
299  f.sw = 0.25 * (x.sw + x.s + x.c + x.w);
300  f.se = 0.25 * (x.s + x.se + x.e + x.c);
301  f.ne = 0.25 * (x.c + x.e + x.ne + x.n);
302  f.nw = 0.25 * (x.w + x.c + x.n + x.nw);
303
304  f.s = 0.5 * (x.c + x.s);
305  f.e = 0.5 * (x.c + x.e);
306  f.n = 0.5 * (x.c + x.n);
307  f.w = 0.5 * (x.c + x.w);
308  }
309
310  // compute the grounding fraction for the current cell by breaking it into 8
311  // triangles
312  double fraction = 0.125 * (grounded_area_fraction(f.c, f.ne, f.n) +
313  grounded_area_fraction(f.c, f.n, f.nw) +
314  grounded_area_fraction(f.c, f.nw, f.w) +
315  grounded_area_fraction(f.c, f.w, f.sw) +
316  grounded_area_fraction(f.c, f.sw, f.s) +
317  grounded_area_fraction(f.c, f.s, f.se) +
318  grounded_area_fraction(f.c, f.se, f.e) +
319  grounded_area_fraction(f.c, f.e, f.ne));
320
321  result(i, j) = clip(fraction, 0.0, 1.0);
322
323  }
324  } catch (...) {
325  loop.failed();
326  }
327  loop.check();
328 }
329
330 } // end of namespace pism
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
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
#define PISM_ERROR_LOCATION
static double F(double SL, double B, double H, double alpha)
Point intersect_bc(double b, double c)
double grounded_area_fraction(double a, double b, double c)
static bool same(const Point &a, const Point &b)
Point intersect_ac(double a, double c)
T clip(T x, T a, T b)
void compute_grounded_cell_fraction(double ice_density, double ocean_density, const array::Scalar1 &sea_level, const array::Scalar1 &ice_thickness, const array::Scalar1 &bed_topography, array::Scalar &result)
stencils::Box< double > Box
bool invalid(const Point &p)
static double triangle_area(const Point &a, const Point &b, const Point &c)
Point intersect_ab(double a, double b)
static double S(unsigned n)
Definition: test_cube.c:58