tgrounded_cell_fraction.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
---
tgrounded_cell_fraction.cc (10069B)
---
1 /* Copyright (C) 2018 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 #include <cmath> // fabs
22
23 #include "grounded_cell_fraction.hh"
24
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/pism_utilities.hh" // clip
27 #include "pism/util/iceModelVec.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 * 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 } else {
75 return {-1.0, -1.0}; // no intersection
76 }
77 }
78
79 /*!
80 * Compute the coordinates of the intersection of (z = 0) with the side BC.
81 */
82 Point intersect_bc(double b, double c) {
83 if (b != c) {
84 return {c / (c - b), b / (b - c)};
85 } else {
86 return {-1.0, -1.0}; // no intersection
87 }
88 }
89
90 /*!
91 * Compute the coordinates of the intersection of (z = 0) with the side AC.
92 */
93 Point intersect_ac(double a, double c) {
94 if (a != c) {
95 return {0.0, a / (a - c)};
96 } else {
97 return {-1.0, -1.0}; // no intersection
98 }
99 }
100
101 /*!
102 * Return true if a point p is not a valid point on a side of the triangle
103 * (0,0)-(1,0)-(0,1).
104 *
105 * This is not a complete test (for example, it does not check if y = 1 - x for points on
106 * the hypotenuse). The point of this is to exclude the kinds of invalid points we are
107 * likely to see, not all of them.
108 *
109 * Note that we use (-1, -1) to indicate "invalid" points in the rest of the code and
110 * these are easy to detect: they require only one comparison.
111 */
112 bool invalid(const Point &p) {
113 if (p.x < 0.0 or p.x > 1.0 or p.y < 0.0 or p.y > 1.0) {
114 return true;
115 } else {
116 return false;
117 }
118 }
119
120 /*!
121 * Return true if two points are the same.
122 */
123 static bool same(const Point &a, const Point &b) {
124 double threshold = 1e-12;
125 return fabs(a.x - b.x) < threshold and fabs(a.y - b.y) < threshold;
126 }
127
128 /*!
129 * Consider the right triangle A(0,0) - B(1,0) - C(0,1).
130 *
131 * Define a linear function z = a + (b - a) * x + (c - a) * y, where a, b, and c are its
132 * values at points A, B, and C, respectively.
133 *
134 * Our goal is to find the fraction of the triangle ABC in where z > 0.
135 *
136 * This corresponds to the grounded area fraction if z is the flotation criterion
137 * function.
138 */
139 double grounded_area_fraction(double a, double b, double c) {
140
141 if (a > 0.0 and b > 0.0 and c > 0.0) {
142 return 1.0;
143 }
144
145 if (a <= 0.0 and b <= 0.0 and c <= 0.0) {
146 return 0.0;
147 }
148
149 // the area of the triangle (0,0)-(1,0)-(0,1)
150 const double total_area = 0.5;
151
152 const Point
153 ab = intersect_ab(a, b),
154 bc = intersect_bc(b, c),
155 ac = intersect_ac(a, c);
156
157 if (invalid(bc)) {
158 assert(not (invalid(ab) or invalid(ac)));
159
160 double ratio = triangle_area({0.0, 0.0}, ab, ac) / total_area;
161 assert((ratio >= 0.0) and (ratio <= 1.0));
162
163 if (a > 0.0) {
164 return ratio;
165 } else {
166 return 1.0 - ratio;
167 }
168 }
169
170 if (invalid(ac)) {
171 assert(not (invalid(ab) or invalid(bc)));
172
173 double ratio = triangle_area({1.0, 0.0}, bc, ab) / total_area;
174 assert((ratio >= 0.0) and (ratio <= 1.0));
175
176 if (b > 0.0) {
177 return ratio;
178 } else {
179 return 1.0 - ratio;
180 }
181 }
182
183 if (invalid(ab)) {
184 assert(not (invalid(bc) or invalid(ac)));
185
186 double ratio = triangle_area({0.0, 1.0}, ac, bc) / total_area;
187 assert((ratio >= 0.0) and (ratio <= 1.0));
188
189 if (c > 0.0) {
190 return ratio;
191 } else {
192 return 1.0 - ratio;
193 }
194 }
195
196 // Note that we know that ab, bc, and ac are all valid.
197
198 // the a == 0 case, the line F = 0 goes through A
199 if (same(ab, ac)) {
200 double ratio = triangle_area({1.0, 0.0}, bc, ab) / total_area;
201 assert((ratio >= 0.0) and (ratio <= 1.0));
202
203 if (b > 0.0) {
204 return ratio;
205 } else {
206 return 1.0 - ratio;
207 }
208 }
209
210 // the b == 0 case and the c == 0 case
211 if (same(ab, bc) or same(ac, bc)) {
212 double ratio = triangle_area({0.0, 0.0}, ab, ac) / total_area;
213 assert((ratio >= 0.0) and (ratio <= 1.0));
214
215 if (a > 0.0) {
216 return ratio;
217 } else {
218 return 1.0 - ratio;
219 }
220 }
221
222 // Note: the case of F=0 coinciding with a side of the triangle is covered by if clauses
223 // above. For example, when F=0 coincides with AC, we have a = c = 0 and intersect_ac(a, c)
224 // returns an invalid intersection point.
225
226 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
227 "the logic in grounded_area_fraction failed! Please submit a bug report.");
228 }
229
230 /*!
231 * The flotation criterion.
232 */
233 static double F(double SL, double B, double H, double alpha) {
234 double
235 water_depth = SL - B,
236 shelf_depth = H * alpha;
237 return shelf_depth - water_depth;
238 }
239
240 typedef BoxStencil<double> Box;
241
242 /*!
243 * Compute the flotation criterion at all the points in the box stencil.
244 */
245 static Box F(const Box &SL, const Box &B, const Box &H, double alpha) {
246 return {F(SL.ij, B.ij, H.ij, alpha),
247 F(SL.n, B.n, H.n, alpha),
248 F(SL.nw, B.nw, H.nw, alpha),
249 F(SL.w, B.w, H.w, alpha),
250 F(SL.sw, B.sw, H.sw, alpha),
251 F(SL.s, B.s, H.s, alpha),
252 F(SL.se, B.se, H.se, alpha),
253 F(SL.e, B.e, H.e, alpha),
254 F(SL.ne, B.ne, H.ne, alpha)};
255 }
256
257 /*!
258 * @param[in] ice_density ice density, kg/m3
259 * @param[in] ocean_density ocean_density, kg/m3
260 * @param[in] sea_level sea level (flotation) elevation, m
261 * @param[in] ice_thickness ice thickness, m
262 * @param[in] bed_topography bed elevation, m
263 * @param[out] result grounded cell fraction, between 0 (floating) and 1 (grounded)
264 */
265 void compute_grounded_cell_fraction(double ice_density,
266 double ocean_density,
267 const IceModelVec2S &sea_level,
268 const IceModelVec2S &ice_thickness,
269 const IceModelVec2S &bed_topography,
270 IceModelVec2S &result) {
271 IceGrid::ConstPtr grid = result.grid();
272 double alpha = ice_density / ocean_density;
273
274 IceModelVec::AccessList list{&sea_level, &ice_thickness, &bed_topography, &result};
275
276 ParallelSection loop(grid->com);
277 try {
278 for (Points p(*grid); p; p.next()) {
279 const int i = p.i(), j = p.j();
280
281 auto S = sea_level.box(i, j);
282 auto H = ice_thickness.box(i, j);
283 auto B = bed_topography.box(i, j);
284
285 auto f = F(S, B, H, alpha);
286
287 /*
288 NW----------------N----------------NE
289 | | |
290 | | |
291 | nw--------n--------ne |
292 | | | | |
293 | | | | |
294 W--------w--------o--------e--------E
295 | | | | |
296 | | | | |
297 | sw--------s--------se |
298 | | |
299 | | |
300 SW----------------S----------------SE
301 */
302
303 double
304 f_o = f.ij,
305 f_sw = 0.25 * (f.sw + f.s + f.ij + f.w),
306 f_se = 0.25 * (f.s + f.se + f.e + f.ij),
307 f_ne = 0.25 * (f.ij + f.e + f.ne + f.n),
308 f_nw = 0.25 * (f.w + f.ij + f.n + f.nw);
309
310 double
311 f_s = 0.5 * (f.ij + f.s),
312 f_e = 0.5 * (f.ij + f.e),
313 f_n = 0.5 * (f.ij + f.n),
314 f_w = 0.5 * (f.ij + f.w);
315
316 double fraction = 0.125 * (grounded_area_fraction(f_o, f_ne, f_n) +
317 grounded_area_fraction(f_o, f_n, f_nw) +
318 grounded_area_fraction(f_o, f_nw, f_w) +
319 grounded_area_fraction(f_o, f_w, f_sw) +
320 grounded_area_fraction(f_o, f_sw, f_s) +
321 grounded_area_fraction(f_o, f_s, f_se) +
322 grounded_area_fraction(f_o, f_se, f_e) +
323 grounded_area_fraction(f_o, f_e, f_ne));
324
325 result(i, j) = clip(fraction, 0.0, 1.0);
326
327 }
328 } catch (...) {
329 loop.failed();
330 }
331 loop.check();
332 }
333
334 } // end of namespace pism