PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
EmptyingProblem.cc
Go to the documentation of this file.
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 
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 
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 /*!
168  * Compute steady state water flux.
169  *
170  * @param[in] geometry ice geometry
171  * @param[in] no_model_mask no model mask
172  * @param[in] water_input_rate water input rate in m/s
173  */
174 void EmptyingProblem::update(const Geometry &geometry,
175  const IceModelVec2Int *no_model_mask,
176  const IceModelVec2S &water_input_rate,
177  bool recompute_potential) {
178 
179  const double
180  eps = 1e-16,
181  cell_area = m_grid->cell_area(),
182  u_max = m_speed,
183  v_max = m_speed,
184  dt = 0.5 / (u_max / m_dx + v_max / m_dy), // CFL condition
185  volume_ratio = m_config->get_number("hydrology.steady.volume_ratio");
186 
187  const int n_iterations = m_config->get_number("hydrology.steady.n_iterations");
188 
189  if (recompute_potential) {
191 
192  compute_mask(geometry.cell_type, no_model_mask, m_domain_mask);
193 
194  // updates ghosts of m_potential
198  m_potential);
199 
200  // diagnostics
201  {
203 
205 
207  }
208  }
209 
210  // set initial state and compute initial volume
211  double volume_0 = 0.0;
212  {
213  IceModelVec::AccessList list{&geometry.cell_type, &m_W, &water_input_rate};
214 
215  for (Points p(*m_grid); p; p.next()) {
216  const int i = p.i(), j = p.j();
217 
218  if (geometry.cell_type.icy(i, j)) {
219  m_W(i, j) = m_tau * water_input_rate(i, j);
220  } else {
221  m_W(i, j) = 0.0;
222  }
223 
224  volume_0 += m_W(i, j);
225  }
226  volume_0 = cell_area * GlobalSum(m_grid->com, volume_0);
227  }
228  m_W.update_ghosts();
229 
230  // uses ghosts of m_potential and m_domain_mask, updates ghosts of m_Vstag
232 
233  m_Qsum.set(0.0);
234 
235  // no input means no flux
236  if (volume_0 == 0.0) {
237  m_Q.set(0.0);
238  m_q_sg.set(0.0);
239  return;
240  }
241 
242  double volume = 0.0;
243  int step_counter = 0;
244 
246 
247  for (step_counter = 0; step_counter < n_iterations; ++step_counter) {
248  volume = 0.0;
249 
250  for (Points p(*m_grid); p; p.next()) {
251  const int i = p.i(), j = p.j();
252 
253  auto v = m_Vstag.star(i, j);
254  auto w = m_W.star(i, j);
255 
256  double
257  q_n = v.n * (v.n >= 0.0 ? w.ij : w.n),
258  q_e = v.e * (v.e >= 0.0 ? w.ij : w.e),
259  q_s = v.s * (v.s >= 0.0 ? w.s : w.ij),
260  q_w = v.w * (v.w >= 0.0 ? w.w : w.ij),
261  divQ = (q_e - q_w) / m_dx + (q_n - q_s) / m_dy;
262 
263  // update water thickness
264  if (m_domain_mask(i, j) > 0.5) {
265  m_tmp(i, j) = w.ij + dt * (- divQ);
266  } else {
267  m_tmp(i, j) = 0.0;
268  }
269 
270  if (m_tmp(i, j) < -eps) {
271  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "W(%d, %d) = %f < 0",
272  i, j, m_tmp(i, j));
273  }
274 
275  // accumulate the water flux
276  m_Qsum(i, j, 0) += dt * q_e;
277  m_Qsum(i, j, 1) += dt * q_n;
278 
279  // compute volume
280  volume += m_tmp(i, j);
281  }
282 
284  volume = cell_area * GlobalSum(m_grid->com, volume);
285 
286  if (volume / volume_0 <= volume_ratio) {
287  break;
288  }
289  m_log->message(3, "%04d V = %f\n", step_counter, volume / volume_0);
290  } // end of the loop
291 
292  m_log->message(3, "Emptying problem: stopped after %d iterations. V = %f\n",
293  step_counter, volume / volume_0);
294 
295  double epsilon = volume / volume_0;
296 
299  true, // include floating ice
300  m_Q);
301  m_Q.scale(1.0 / (m_tau * (1.0 - epsilon)));
302 
304 }
305 
306 /*! Compute the unmodified hydraulic potential (with sinks).
307  *
308  * @param[in] H ice thickness
309  * @param[in] b ice bottom surface elevation
310  * @param[out] result simplified hydraulic potential used by to compute velocity
311  */
313  const IceModelVec2S &b,
314  IceModelVec2S &result) const {
315  const double
316  g = m_config->get_number("constants.standard_gravity"),
317  rho_i = m_config->get_number("constants.ice.density"),
318  rho_w = m_config->get_number("constants.fresh_water.density");
319 
320  IceModelVec::AccessList list({&H, &b, &result});
321 
322  for (Points p(*m_grid); p; p.next()) {
323  const int i = p.i(), j = p.j();
324 
325  result(i, j) = rho_i * g * H(i, j) + rho_w * g * b(i, j);
326  }
327 
328  result.update_ghosts();
329 }
330 
333  const IceModelVec2Int &domain_mask,
334  IceModelVec2S &result) {
335  IceModelVec2S &psi_new = m_tmp;
336 
337  double delta = m_config->get_number("hydrology.steady.potential_delta");
338 
339  int n_iterations = m_config->get_number("hydrology.steady.potential_n_iterations");
340  int step_counter = 0;
341  int n_sinks = 0;
342  int n_sinks_remaining = 0;
343 
344  compute_raw_potential(ice_thickness, ice_bottom_surface, result);
345 
346  IceModelVec::AccessList list{&result, &psi_new, &domain_mask};
347  for (step_counter = 0; step_counter < n_iterations; ++step_counter) {
348 
349  n_sinks_remaining = 0;
350  for (Points p(*m_grid); p; p.next()) {
351  const int i = p.i(), j = p.j();
352 
353  if (domain_mask(i, j) > 0.5) {
354  auto P = result.star(i, j);
355 
356  double
357  v_n = - (P.n - P.ij),
358  v_e = - (P.e - P.ij),
359  v_s = - (P.ij - P.s),
360  v_w = - (P.ij - P.w);
361 
362  if (v_e <= 0.0 and v_w >= 0.0 and v_n <= 0.0 and v_s >= 0.0) {
363  ++n_sinks_remaining;
364 
365  psi_new(i, j) = 0.25 * (P.n + P.e + P.s + P.w) + delta;
366  } else {
367  psi_new(i, j) = result(i, j);
368  }
369  } else {
370  psi_new(i, j) = result(i, j);
371  }
372  }
373  // this call updates ghosts of result
374  result.copy_from(psi_new);
375 
376  n_sinks_remaining = GlobalSum(m_grid->com, n_sinks_remaining);
377 
378  // remember the original number of sinks
379  if (step_counter == 0) {
380  n_sinks = n_sinks_remaining;
381  }
382 
383  if (n_sinks_remaining == 0) {
384  m_log->message(3, "Emptying problem: filled %d sinks after %d iterations.\n",
385  n_sinks - n_sinks_remaining, step_counter);
386  break;
387  }
388  } // end of loop over step_counter
389 
390  if (n_sinks_remaining > 0) {
391  m_log->message(2, "WARNING: %d sinks left.\n", n_sinks_remaining);
392  }
393 }
394 
395 
396 static double K(double psi_x, double psi_y, double speed, double epsilon) {
397  return speed / std::max(Vector2(psi_x, psi_y).magnitude(), epsilon);
398 }
399 
400 /*!
401  * Compute water velocity on the staggered grid.
402  */
404  const IceModelVec2Int &domain_mask,
405  IceModelVec2Stag &result) const {
406 
407  IceModelVec::AccessList list{&psi, &result, &domain_mask};
408 
409  for (Points p(*m_grid); p; p.next()) {
410  const int i = p.i(), j = p.j();
411 
412  for (int o = 0; o < 2; ++o) {
413  double psi_x, psi_y;
414  if (o == 0) {
415  psi_x = (psi(i + 1, j) - psi(i, j)) / m_dx;
416  psi_y = (psi(i + 1, j + 1) + psi(i, j + 1) - psi(i + 1, j - 1) - psi(i, j - 1)) / (4.0 * m_dy);
417  result(i, j, o) = - K(psi_x, psi_y, m_speed, m_eps_gradient) * psi_x;
418  } else {
419  psi_x = (psi(i + 1, j + 1) + psi(i + 1, j) - psi(i - 1, j + 1) - psi(i - 1, j)) / (4.0 * m_dx);
420  psi_y = (psi(i, j + 1) - psi(i, j)) / m_dy;
421  result(i, j, o) = - K(psi_x, psi_y, m_speed, m_eps_gradient) * psi_y;
422  }
423 
424  auto M = domain_mask.star(i, j);
425 
426  if (M.ij == 0 and M.e == 0) {
427  result(i, j, 0) = 0.0;
428  }
429 
430  if (M.ij == 0 and M.n == 0) {
431  result(i, j, 1) = 0.0;
432  }
433  }
434  }
435  result.update_ghosts();
436 }
437 
438 /*!
439  * Compute the mask that defines the domain: ones in the domain, zeroes elsewhere.
440  */
442  const IceModelVec2Int *no_model_mask,
443  IceModelVec2Int &result) const {
444 
445  IceModelVec::AccessList list{&cell_type, &result};
446 
447  if (no_model_mask) {
448  list.add(*no_model_mask);
449  }
450 
451  for (Points p(*m_grid); p; p.next()) {
452  const int i = p.i(), j = p.j();
453 
454  if (not cell_type.ice_free_ocean(i, j)) {
455  result(i, j) = 1.0;
456  } else {
457  result(i, j) = 0.0;
458  }
459 
460  if (no_model_mask and no_model_mask->as_int(i, j) == 1) {
461  result(i, j) = 0.0;
462  }
463  }
464 
465  result.update_ghosts();
466 }
467 
468 
469 
471  return {{"steady_state_hydraulic_potential", Diagnostic::wrap(m_potential)},
472  {"hydraulic_potential_adjustment", Diagnostic::wrap(m_adjustment)},
473  {"hydraulic_sinks", Diagnostic::wrap(m_sinks)},
474  {"remaining_water_thickness", Diagnostic::wrap(m_W)},
475  {"effective_water_velocity", Diagnostic::wrap(m_q_sg)}};
476 }
477 
478 
479 /*!
480  * Remaining water thickness.
481  *
482  * This field can be used to get an idea of where water is accumulated. (This affects the
483  * quality of the estimate of the water flux).
484  */
486  return m_W;
487 }
488 
489 /*!
490  * Steady state water flux.
491  */
493  return m_Q;
494 }
495 
496 /*!
497  * Effective water velocity (flux per unit area of the front).
498  */
500  return m_q_sg;
501 }
502 
503 /*!
504  * Hydraulic potential used to determine flow direction.
505  */
507  return m_potential;
508 }
509 
510 /*!
511  * Map of sinks.
512  */
514  return m_sinks;
515 }
516 
517 /*!
518  * Adjustment applied to the unmodified hydraulic potential to eliminate sinks.
519  */
521  return m_adjustment;
522 }
523 
524 } // end of namespace hydrology
525 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:101
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S sea_level_elevation
Definition: Geometry.hh:50
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool icy(int i, int j) const
bool ice_free_ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
int as_int(int i, int j) const
stencils::Star< int > star(int i, int j) const
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
void add(double alpha, const IceModelVec2S &x)
void copy_from(const IceModelVec2S &source)
stencils::Star< double > star(int i, int j) const
stencils::Star< double > star(int i, int j) const
Returns the values at interfaces of the cell i,j using the staggered grid.
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
const IceModelVec2V & flux() const
virtual void compute_raw_potential(const IceModelVec2S &ice_thickness, const IceModelVec2S &ice_bottom_surface, IceModelVec2S &result) const
void update(const Geometry &geometry, const IceModelVec2Int *no_model_mask, const IceModelVec2S &water_input_rate, bool recompute_potential=true)
const IceModelVec2Int & sinks() const
DiagnosticList diagnostics() const
EmptyingProblem(IceGrid::ConstPtr g)
const IceModelVec2V & effective_water_velocity() const
const IceModelVec2S & potential() const
void compute_velocity(const IceModelVec2S &hydraulic_potential, const IceModelVec2Int &mask, IceModelVec2Stag &result) const
void compute_potential(const IceModelVec2S &ice_thickness, const IceModelVec2S &ice_bottom_surface, const IceModelVec2Int &domain_mask, IceModelVec2S &result)
const IceModelVec2S & remaining_water_thickness() const
void compute_mask(const IceModelVec2CellType &cell_type, const IceModelVec2Int *no_model_mask, IceModelVec2Int &result) const
const IceModelVec2S & adjustment() const
#define PISM_ERROR_LOCATION
static void effective_water_velocity(const Geometry &geometry, const IceModelVec2V &water_flux, IceModelVec2V &result)
static void compute_sinks(const IceModelVec2Int &domain_mask, const IceModelVec2S &psi, IceModelVec2S &result)
static double K(double psi_x, double psi_y, double speed, double epsilon)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double g
Definition: exactTestP.cc:39
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result)
Definition: Geometry.cc:223
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
void staggered_to_regular(const IceModelVec2CellType &cell_type, const IceModelVec2Stag &input, bool include_floating_ice, IceModelVec2S &result)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49