tEmptyingProblem.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
---
tEmptyingProblem.cc (16490B)
---
1 /* Copyright (C) 2019, 2020 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 "EmptyingProblem.hh"
21
22 #include "pism/geometry/Geometry.hh"
23 #include "pism/util/pism_utilities.hh"
24 #include "pism/geometry/Geometry.hh"
25
26 namespace pism {
27 namespace hydrology {
28
29 namespace diagnostics {
30
31 /*! Compute the map of sinks.
32 *
33 * This field is purely diagnostic: it can be used to diagnose failures of the code
34 * filling dips to eliminate sinks (and get a better estimate of the steady state flow
35 * direction).
36 */
37 static void compute_sinks(const IceModelVec2Int &domain_mask,
38 const IceModelVec2S &psi,
39 IceModelVec2S &result) {
40
41 IceGrid::ConstPtr grid = result.grid();
42
43 IceModelVec::AccessList list{&psi, &domain_mask, &result};
44
45 for (Points p(*grid); p; p.next()) {
46 const int i = p.i(), j = p.j();
47
48 auto P = psi.star(i, j);
49
50 double
51 v_n = - (P.n - P.ij),
52 v_e = - (P.e - P.ij),
53 v_s = - (P.ij - P.s),
54 v_w = - (P.ij - P.w);
55
56 if (domain_mask(i, j) > 0.5 and v_e <= 0.0 and v_w >= 0.0 and v_n <= 0.0 and v_s >= 0.0) {
57 result(i, j) = 1.0;
58 } else {
59 result(i, j) = 0.0;
60 }
61 }
62 }
63
64 static void effective_water_velocity(const Geometry &geometry,
65 const IceModelVec2V &water_flux,
66 IceModelVec2V &result) {
67
68 IceGrid::ConstPtr grid = result.grid();
69
70 const IceModelVec2CellType &cell_type = geometry.cell_type;
71 const IceModelVec2S &bed_elevation = geometry.bed_elevation;
72 const IceModelVec2S &ice_thickness = geometry.ice_thickness;
73 const IceModelVec2S &sea_level_elevation = geometry.sea_level_elevation;
74
75 IceModelVec::AccessList list
76 {&ice_thickness, &bed_elevation, &cell_type, &sea_level_elevation,
77 &water_flux, &result};
78
79 double
80 grid_spacing = 0.5 * (grid->dx() + grid->dy()),
81 eps = 1.0; // q_sg regularization
82
83 for (Points p(*grid); p; p.next()) {
84 const int i = p.i(), j = p.j();
85
86 if (cell_type.icy(i, j)) {
87 // Convert subglacial water flux (m^2/s) to an "effective subglacial freshwater
88 // velocity" or flux per unit area of ice front in m/day (see Xu et al 2013, section
89 // 2, paragraph 11).
90 //
91 // [flux] = m^2 / s, so
92 // [flux * grid_spacing] = m^3 / s, so
93 // [flux * grid_spacing / submerged_front_area] = m / s, and
94 // [flux * grid_spacing * (s / day) / submerged_front_area] = m / day
95
96 double water_depth = std::max(sea_level_elevation(i, j) - bed_elevation(i, j), 0.0),
97 submerged_front_area = water_depth * grid_spacing;
98
99 if (submerged_front_area > 0.0) {
100 auto Q_sg = water_flux(i, j) * grid_spacing;
101 auto q_sg = Q_sg / std::max(submerged_front_area, eps);
102
103 result(i, j) = q_sg;
104 } else {
105 result(i, j) = 0.0;
106 }
107 } else {
108 result(i, j) = 0.0;
109 }
110 } // end of the loop over grid points
111 }
112
113 } // end of namespace diagnostics
114
115 EmptyingProblem::EmptyingProblem(IceGrid::ConstPtr grid)
116 : Component(grid),
117 m_potential(grid, "hydraulic_potential", WITH_GHOSTS, 1),
118 m_tmp(grid, "temporary_storage", WITHOUT_GHOSTS),
119 m_bottom_surface(grid, "ice_bottom_surface", WITHOUT_GHOSTS),
120 m_W(grid, "remaining_water_thickness", WITH_GHOSTS, 1),
121 m_Vstag(grid, "V_staggered", WITH_GHOSTS),
122 m_Qsum(grid, "flux_total", WITH_GHOSTS, 1),
123 m_domain_mask(grid, "domain_mask", WITH_GHOSTS, 1),
124 m_Q(grid, "_water_flux", WITHOUT_GHOSTS),
125 m_q_sg(grid, "_effective_water_velocity", WITHOUT_GHOSTS),
126 m_adjustment(grid, "hydraulic_potential_adjustment", WITHOUT_GHOSTS),
127 m_sinks(grid, "sinks", WITHOUT_GHOSTS) {
128
129 m_potential.set_attrs("diagnostic", "estimate of the steady state hydraulic potential in the steady hydrology model",
130 "Pa", "Pa", "", 0);
131
132 m_bottom_surface.set_attrs("internal", "ice bottom surface elevation",
133 "m", "m", "", 0);
134
135 m_W.set_attrs("diagnostic",
136 "scaled water thickness in the steady state hydrology model"
137 " (has no physical meaning)",
138 "m", "m", "", 0);
139
140 m_Vstag.set_attrs("diagnostic", "water velocity on the staggered grid",
141 "m s-1", "m s-1", "", 0);
142
143 m_domain_mask.set_attrs("internal", "mask defining the domain", "", "", "", 0);
144
145 m_Q.set_attrs("diagnostic", "steady state water flux", "m2 s-1", "m2 s-1", "", 0);
146
147 m_q_sg.set_attrs("diagnostic", "x-component of the effective water velocity in the steady-state hydrology model",
148 "m s-1", "m day-1", "", 0);
149 m_q_sg.set_attrs("diagnostic", "y-component of the effective water velocity in the steady-state hydrology model",
150 "m s-1", "m day-1", "", 1);
151
152 m_sinks.set_attrs("diagnostic", "map of sinks in the domain (for debugging)", "", "", "", 0);
153
154 m_adjustment.set_attrs("diagnostic",
155 "potential adjustment needed to fill sinks"
156 " when computing an estimate of the steady-state hydraulic potential",
157 "Pa", "Pa", "", 0);
158
159 m_eps_gradient = 1e-2;
160 m_speed = 1.0;
161
162 m_dx = m_grid->dx();
163 m_dy = m_grid->dy();
164 m_tau = m_config->get_number("hydrology.steady.input_rate_scaling");
165 }
166
167 EmptyingProblem::~EmptyingProblem() {
168 // empty
169 }
170
171 /*!
172 * Compute steady state water flux.
173 *
174 * @param[in] geometry ice geometry
175 * @param[in] no_model_mask no model mask
176 * @param[in] water_input_rate water input rate in m/s
177 */
178 void EmptyingProblem::update(const Geometry &geometry,
179 const IceModelVec2Int *no_model_mask,
180 const IceModelVec2S &water_input_rate,
181 bool recompute_potential) {
182
183 const double
184 eps = 1e-16,
185 cell_area = m_grid->cell_area(),
186 u_max = m_speed,
187 v_max = m_speed,
188 dt = 0.5 / (u_max / m_dx + v_max / m_dy), // CFL condition
189 volume_ratio = m_config->get_number("hydrology.steady.volume_ratio");
190
191 const int n_iterations = m_config->get_number("hydrology.steady.n_iterations");
192
193 if (recompute_potential) {
194 ice_bottom_surface(geometry, m_bottom_surface);
195
196 compute_mask(geometry.cell_type, no_model_mask, m_domain_mask);
197
198 // updates ghosts of m_potential
199 compute_potential(geometry.ice_thickness,
200 m_bottom_surface,
201 m_domain_mask,
202 m_potential);
203
204 // diagnostics
205 {
206 compute_raw_potential(geometry.ice_thickness, m_bottom_surface, m_adjustment);
207
208 m_potential.add(-1.0, m_adjustment, m_adjustment);
209
210 diagnostics::compute_sinks(m_domain_mask, m_potential, m_sinks);
211 }
212 }
213
214 // set initial state and compute initial volume
215 double volume_0 = 0.0;
216 {
217 IceModelVec::AccessList list{&geometry.cell_type, &m_W, &water_input_rate};
218
219 for (Points p(*m_grid); p; p.next()) {
220 const int i = p.i(), j = p.j();
221
222 if (geometry.cell_type.icy(i, j)) {
223 m_W(i, j) = m_tau * water_input_rate(i, j);
224 } else {
225 m_W(i, j) = 0.0;
226 }
227
228 volume_0 += m_W(i, j);
229 }
230 volume_0 = cell_area * GlobalSum(m_grid->com, volume_0);
231 }
232 m_W.update_ghosts();
233
234 // uses ghosts of m_potential and m_domain_mask, updates ghosts of m_Vstag
235 compute_velocity(m_potential, m_domain_mask, m_Vstag);
236
237 m_Qsum.set(0.0);
238
239 // no input means no flux
240 if (volume_0 == 0.0) {
241 m_Q.set(0.0);
242 m_q_sg.set(0.0);
243 return;
244 }
245
246 double volume = 0.0;
247 int step_counter = 0;
248
249 IceModelVec::AccessList list{&m_Qsum, &m_W, &m_Vstag, &m_domain_mask, &m_tmp};
250
251 for (step_counter = 0; step_counter < n_iterations; ++step_counter) {
252 volume = 0.0;
253
254 for (Points p(*m_grid); p; p.next()) {
255 const int i = p.i(), j = p.j();
256
257 auto v = m_Vstag.star(i, j);
258 auto w = m_W.star(i, j);
259
260 double
261 q_n = v.n * (v.n >= 0.0 ? w.ij : w.n),
262 q_e = v.e * (v.e >= 0.0 ? w.ij : w.e),
263 q_s = v.s * (v.s >= 0.0 ? w.s : w.ij),
264 q_w = v.w * (v.w >= 0.0 ? w.w : w.ij),
265 divQ = (q_e - q_w) / m_dx + (q_n - q_s) / m_dy;
266
267 // update water thickness
268 if (m_domain_mask(i, j) > 0.5) {
269 m_tmp(i, j) = w.ij + dt * (- divQ);
270 } else {
271 m_tmp(i, j) = 0.0;
272 }
273
274 if (m_tmp(i, j) < -eps) {
275 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "W(%d, %d) = %f < 0",
276 i, j, m_tmp(i, j));
277 }
278
279 // accumulate the water flux
280 m_Qsum(i, j, 0) += dt * q_e;
281 m_Qsum(i, j, 1) += dt * q_n;
282
283 // compute volume
284 volume += m_tmp(i, j);
285 }
286
287 m_W.copy_from(m_tmp);
288 volume = cell_area * GlobalSum(m_grid->com, volume);
289
290 if (volume / volume_0 <= volume_ratio) {
291 break;
292 }
293 m_log->message(3, "%04d V = %f\n", step_counter, volume / volume_0);
294 } // end of the loop
295
296 m_log->message(3, "Emptying problem: stopped after %d iterations. V = %f\n",
297 step_counter, volume / volume_0);
298
299 double epsilon = volume / volume_0;
300
301 m_Qsum.update_ghosts();
302 staggered_to_regular(geometry.cell_type, m_Qsum,
303 true, // include floating ice
304 m_Q);
305 m_Q.scale(1.0 / (m_tau * (1.0 - epsilon)));
306
307 diagnostics::effective_water_velocity(geometry, m_Q, m_q_sg);
308 }
309
310 /*! Compute the unmodified hydraulic potential (with sinks).
311 *
312 * @param[in] H ice thickness
313 * @param[in] b ice bottom surface elevation
314 * @param[out] result simplified hydraulic potential used by to compute velocity
315 */
316 void EmptyingProblem::compute_raw_potential(const IceModelVec2S &H,
317 const IceModelVec2S &b,
318 IceModelVec2S &result) const {
319 const double
320 g = m_config->get_number("constants.standard_gravity"),
321 rho_i = m_config->get_number("constants.ice.density"),
322 rho_w = m_config->get_number("constants.fresh_water.density");
323
324 IceModelVec::AccessList list({&H, &b, &result});
325
326 for (Points p(*m_grid); p; p.next()) {
327 const int i = p.i(), j = p.j();
328
329 result(i, j) = rho_i * g * H(i, j) + rho_w * g * b(i, j);
330 }
331
332 result.update_ghosts();
333 }
334
335 void EmptyingProblem::compute_potential(const IceModelVec2S &ice_thickness,
336 const IceModelVec2S &ice_bottom_surface,
337 const IceModelVec2Int &domain_mask,
338 IceModelVec2S &result) {
339 IceModelVec2S &psi_new = m_tmp;
340
341 double delta = m_config->get_number("hydrology.steady.potential_delta");
342
343 int n_iterations = m_config->get_number("hydrology.steady.potential_n_iterations");
344 int step_counter = 0;
345 int n_sinks = 0;
346 int n_sinks_remaining = 0;
347
348 compute_raw_potential(ice_thickness, ice_bottom_surface, result);
349
350 IceModelVec::AccessList list{&result, &psi_new, &domain_mask};
351 for (step_counter = 0; step_counter < n_iterations; ++step_counter) {
352
353 n_sinks_remaining = 0;
354 for (Points p(*m_grid); p; p.next()) {
355 const int i = p.i(), j = p.j();
356
357 if (domain_mask(i, j) > 0.5) {
358 auto P = result.star(i, j);
359
360 double
361 v_n = - (P.n - P.ij),
362 v_e = - (P.e - P.ij),
363 v_s = - (P.ij - P.s),
364 v_w = - (P.ij - P.w);
365
366 if (v_e <= 0.0 and v_w >= 0.0 and v_n <= 0.0 and v_s >= 0.0) {
367 ++n_sinks_remaining;
368
369 psi_new(i, j) = 0.25 * (P.n + P.e + P.s + P.w) + delta;
370 } else {
371 psi_new(i, j) = result(i, j);
372 }
373 } else {
374 psi_new(i, j) = result(i, j);
375 }
376 }
377 // this call updates ghosts of result
378 result.copy_from(psi_new);
379
380 n_sinks_remaining = GlobalSum(m_grid->com, n_sinks_remaining);
381
382 // remember the original number of sinks
383 if (step_counter == 0) {
384 n_sinks = n_sinks_remaining;
385 }
386
387 if (n_sinks_remaining == 0) {
388 m_log->message(3, "Emptying problem: filled %d sinks after %d iterations.\n",
389 n_sinks - n_sinks_remaining, step_counter);
390 break;
391 }
392 } // end of loop over step_counter
393
394 if (n_sinks_remaining > 0) {
395 m_log->message(2, "WARNING: %d sinks left.\n", n_sinks_remaining);
396 }
397 }
398
399
400 static double K(double psi_x, double psi_y, double speed, double epsilon) {
401 return speed / std::max(Vector2(psi_x, psi_y).magnitude(), epsilon);
402 }
403
404 /*!
405 * Compute water velocity on the staggered grid.
406 */
407 void EmptyingProblem::compute_velocity(const IceModelVec2S &psi,
408 const IceModelVec2Int &domain_mask,
409 IceModelVec2Stag &result) const {
410
411 IceModelVec::AccessList list{&psi, &result, &domain_mask};
412
413 for (Points p(*m_grid); p; p.next()) {
414 const int i = p.i(), j = p.j();
415
416 for (int o = 0; o < 2; ++o) {
417 double psi_x, psi_y;
418 if (o == 0) {
419 psi_x = (psi(i + 1, j) - psi(i, j)) / m_dx;
420 psi_y = (psi(i + 1, j + 1) + psi(i, j + 1) - psi(i + 1, j - 1) - psi(i, j - 1)) / (4.0 * m_dy);
421 result(i, j, o) = - K(psi_x, psi_y, m_speed, m_eps_gradient) * psi_x;
422 } else {
423 psi_x = (psi(i + 1, j + 1) + psi(i + 1, j) - psi(i - 1, j + 1) - psi(i - 1, j)) / (4.0 * m_dx);
424 psi_y = (psi(i, j + 1) - psi(i, j)) / m_dy;
425 result(i, j, o) = - K(psi_x, psi_y, m_speed, m_eps_gradient) * psi_y;
426 }
427
428 auto M = domain_mask.int_star(i, j);
429
430 if (M.ij == 0 and M.e == 0) {
431 result(i, j, 0) = 0.0;
432 }
433
434 if (M.ij == 0 and M.n == 0) {
435 result(i, j, 1) = 0.0;
436 }
437 }
438 }
439 result.update_ghosts();
440 }
441
442 /*!
443 * Compute the mask that defines the domain: ones in the domain, zeroes elsewhere.
444 */
445 void EmptyingProblem::compute_mask(const IceModelVec2CellType &cell_type,
446 const IceModelVec2Int *no_model_mask,
447 IceModelVec2Int &result) const {
448
449 IceModelVec::AccessList list{&cell_type, &result};
450
451 if (no_model_mask) {
452 list.add(*no_model_mask);
453 }
454
455 for (Points p(*m_grid); p; p.next()) {
456 const int i = p.i(), j = p.j();
457
458 if (not cell_type.ice_free_ocean(i, j)) {
459 result(i, j) = 1.0;
460 } else {
461 result(i, j) = 0.0;
462 }
463
464 if (no_model_mask and no_model_mask->as_int(i, j) == 1) {
465 result(i, j) = 0.0;
466 }
467 }
468
469 result.update_ghosts();
470 }
471
472
473
474 DiagnosticList EmptyingProblem::diagnostics() const {
475 return {{"steady_state_hydraulic_potential", Diagnostic::wrap(m_potential)},
476 {"hydraulic_potential_adjustment", Diagnostic::wrap(m_adjustment)},
477 {"hydraulic_sinks", Diagnostic::wrap(m_sinks)},
478 {"remaining_water_thickness", Diagnostic::wrap(m_W)},
479 {"effective_water_velocity", Diagnostic::wrap(m_q_sg)}};
480 }
481
482
483 /*!
484 * Remaining water thickness.
485 *
486 * This field can be used to get an idea of where water is accumulated. (This affects the
487 * quality of the estimate of the water flux).
488 */
489 const IceModelVec2S& EmptyingProblem::remaining_water_thickness() const {
490 return m_W;
491 }
492
493 /*!
494 * Steady state water flux.
495 */
496 const IceModelVec2V& EmptyingProblem::flux() const {
497 return m_Q;
498 }
499
500 /*!
501 * Effective water velocity (flux per unit area of the front).
502 */
503 const IceModelVec2V& EmptyingProblem::effective_water_velocity() const {
504 return m_q_sg;
505 }
506
507 /*!
508 * Hydraulic potential used to determine flow direction.
509 */
510 const IceModelVec2S& EmptyingProblem::potential() const {
511 return m_potential;
512 }
513
514 /*!
515 * Map of sinks.
516 */
517 const IceModelVec2Int& EmptyingProblem::sinks() const {
518 return m_sinks;
519 }
520
521 /*!
522 * Adjustment applied to the unmodified hydraulic potential to eliminate sinks.
523 */
524 const IceModelVec2S& EmptyingProblem::adjustment() const {
525 return m_adjustment;
526 }
527
528 } // end of namespace hydrology
529 } // end of namespace pism