PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
SIAFD.cc
Go to the documentation of this file.
1 // Copyright (C) 2004--2021 Jed Brown, Craig Lingle, Ed Bueler and Constantine Khroulev
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 #include <cstdlib>
20 #include <cassert>
21 
22 #include "SIAFD.hh"
23 #include "BedSmoother.hh"
24 #include "pism/util/EnthalpyConverter.hh"
25 #include "pism/rheology/FlowLawFactory.hh"
26 #include "pism/rheology/grain_size_vostok.hh"
27 #include "pism/util/IceGrid.hh"
28 #include "pism/util/Mask.hh"
29 #include "pism/util/Vars.hh"
30 #include "pism/util/error_handling.hh"
31 #include "pism/util/pism_utilities.hh"
32 #include "pism/util/Profiling.hh"
33 #include "pism/util/IceModelVec2CellType.hh"
34 #include "pism/geometry/Geometry.hh"
35 #include "pism/stressbalance/StressBalance.hh"
36 
37 #include "pism/util/Time.hh"
38 #include "pism/util/Context.hh"
39 
40 namespace pism {
41 namespace stressbalance {
42 
44  : SSB_Modifier(std::move(g)),
45  m_stencil_width(m_config->get_number("grid.max_stencil_width")),
46  m_work_2d_0(m_grid, "work_vector_2d_0", WITH_GHOSTS, m_stencil_width),
47  m_work_2d_1(m_grid, "work_vector_2d_1", WITH_GHOSTS, m_stencil_width),
48  m_h_x(m_grid, "h_x", WITH_GHOSTS),
49  m_h_y(m_grid, "h_y", WITH_GHOSTS),
50  m_D(m_grid, "diffusivity", WITH_GHOSTS),
51  m_delta_0(m_grid, "delta_0", WITH_GHOSTS, m_grid->z()),
52  m_delta_1(m_grid, "delta_1", WITH_GHOSTS, m_grid->z()),
53  m_work_3d_0(m_grid, "work_3d_0", WITH_GHOSTS, m_grid->z()),
54  m_work_3d_1(m_grid, "work_3d_1", WITH_GHOSTS, m_grid->z())
55 {
56  // bed smoother
58 
59  m_seconds_per_year = units::convert(m_sys, 1, "second", "years");
60 
61  m_e_factor = m_config->get_number("stress_balance.sia.enhancement_factor");
62  m_e_factor_interglacial = m_config->get_number("stress_balance.sia.enhancement_factor_interglacial");
63 
64  {
65  rheology::FlowLawFactory ice_factory("stress_balance.sia.", m_config, m_EC);
66  m_flow_law = ice_factory.create();
67  }
68 
69  const bool compute_grain_size_using_age = m_config->get_flag("stress_balance.sia.grain_size_age_coupling");
70  const bool age_model_enabled = m_config->get_flag("age.enabled");
71  const bool e_age_coupling = m_config->get_flag("stress_balance.sia.e_age_coupling");
72 
73  if (compute_grain_size_using_age) {
74  if (not FlowLawUsesGrainSize(*m_flow_law)) {
75  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "flow law %s does not use grain size "
76  "but sia.grain_size_age_coupling was set",
77  m_flow_law->name().c_str());
78  }
79 
80  if (not age_model_enabled) {
81  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "SIAFD: age model is not active but\n"
82  "age is needed for grain-size-based flow law %s",
83  m_flow_law->name().c_str());
84  }
85  }
86 
87  if (e_age_coupling and not age_model_enabled) {
88  throw RuntimeError(PISM_ERROR_LOCATION, "SIAFD: age model is not active but\n"
89  "age is needed for age-dependent flow enhancement");
90  }
91 
92  m_eemian_start = m_config->get_number("time.eemian_start", "seconds");
93  m_eemian_end = m_config->get_number("time.eemian_end", "seconds");
94  m_holocene_start = m_config->get_number("time.holocene_start", "seconds");
95 }
96 
98  delete m_bed_smoother;
99 }
100 
101 //! \brief Initialize the SIA module.
102 void SIAFD::init() {
103 
105 
106  m_log->message(2,
107  "* Initializing the SIA stress balance modifier...\n");
108  m_log->message(2,
109  " [using the %s flow law]\n", m_flow_law->name().c_str());
110 
111 
112  // implements an option e.g. described in @ref Greve97Greenland that is the
113  // enhancement factor is coupled to the age of the ice
114  if (m_config->get_flag("stress_balance.sia.e_age_coupling")) {
115  m_log->message(2,
116  " using age-dependent enhancement factor:\n"
117  " e=%f for ice accumulated during interglacial periods\n"
118  " e=%f for ice accumulated during glacial periods\n",
120  }
121 }
122 
123 //! \brief Do the update; if full_update == false skip the update of 3D velocities and strain
124 //! heating.
125 void SIAFD::update(const IceModelVec2V &sliding_velocity,
126  const Inputs &inputs,
127  bool full_update) {
128 
129  const Profiling &profiling = m_grid->ctx()->profiling();
130 
131  // Check if the smoothed bed computed by BedSmoother is out of date and
132  // recompute if necessary.
133  if (inputs.new_bed_elevation) {
134  profiling.begin("sia.bed_smoother");
136  profiling.end("sia.bed_smoother");
137  }
138 
139  profiling.begin("sia.gradient");
141  profiling.end("sia.gradient");
142 
143  profiling.begin("sia.flux");
144  compute_diffusivity(full_update,
145  *inputs.geometry,
146  inputs.enthalpy,
147  inputs.age,
148  m_h_x, m_h_y, m_D);
150  profiling.end("sia.flux");
151 
152  if (full_update) {
153  profiling.begin("sia.3d_velocity");
154  compute_3d_horizontal_velocity(*inputs.geometry, m_h_x, m_h_y, sliding_velocity,
155  m_u, m_v);
156  profiling.end("sia.3d_velocity");
157  }
158 }
159 
160 
161 //! \brief Compute the ice surface gradient for the SIA.
162 /*!
163  There are three methods for computing the surface gradient. Which method is
164  controlled by configuration parameter `sia.surface_gradient_method` which can
165  have values `haseloff`, `mahaffy`, or `eta`.
166 
167  The most traditional method is to directly differentiate the surface
168  elevation \f$h\f$ by the Mahaffy method [\ref Mahaffy]. The `haseloff` method,
169  suggested by Marianne Haseloff, modifies the Mahaffy method only where
170  ice-free adjacent bedrock points are above the ice surface, and in those
171  cases the returned gradient component is zero.
172 
173  The alternative method, when `sia.surface_gradient_method` = `eta`, transforms
174  the thickness to something more regular and differentiates that. We get back
175  to the gradient of the surface by applying the chain rule. In particular, as
176  shown in [\ref Vazquezetal2003] for the flat bed and \f$n=3\f$ case, if we define
177 
178  \f[\eta = H^{(2n+2)/n}\f]
179 
180  then \f$\eta\f$ is more regular near the margin than \f$H\f$. So we compute
181  the surface gradient by
182 
183  \f[\nabla h = \frac{n}{(2n+2)} \eta^{(-n-2)/(2n+2)} \nabla \eta + \nabla b,\f]
184 
185  recalling that \f$h = H + b\f$. This method is only applied when \f$\eta >
186  0\f$ at a given point; otherwise \f$\nabla h = \nabla b\f$.
187 
188  In all cases we are computing the gradient by finite differences onto a
189  staggered grid. In the method with \f$\eta\f$ we apply centered differences
190  using (roughly) the same method for \f$\eta\f$ and \f$b\f$ that applies
191  directly to the surface elevation \f$h\f$ in the `mahaffy` and `haseloff`
192  methods.
193 
194  \param[out] h_x the X-component of the surface gradient, on the staggered grid
195  \param[out] h_y the Y-component of the surface gradient, on the staggered grid
196 */
198  IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) {
199 
200  const std::string method = m_config->get_string("stress_balance.sia.surface_gradient_method");
201 
202  if (method == "eta") {
203 
205  inputs.geometry->bed_elevation,
206  h_x, h_y);
207 
208  } else if (method == "haseloff") {
209 
211  inputs.geometry->cell_type,
212  h_x, h_y);
213 
214  } else if (method == "mahaffy") {
215 
217  h_x, h_y);
218 
219  } else {
221  "value of sia.surface_gradient_method, option '-gradient %s', is not valid",
222  method.c_str());
223  }
224 }
225 
226 //! \brief Compute the ice surface gradient using the eta-transformation.
227 void SIAFD::surface_gradient_eta(const IceModelVec2S &ice_thickness,
228  const IceModelVec2S &bed_elevation,
229  IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) {
230  const double n = m_flow_law->exponent(), // presumably 3.0
231  etapow = (2.0 * n + 2.0)/n, // = 8/3 if n = 3
232  invpow = 1.0 / etapow,
233  dinvpow = (- n - 2.0) / (2.0 * n + 2.0);
234  const double dx = m_grid->dx(), dy = m_grid->dy(); // convenience
235  IceModelVec2S &eta = m_work_2d_0;
236 
237  // compute eta = H^{8/3}, which is more regular, on reg grid
238 
239  IceModelVec::AccessList list{&eta, &ice_thickness, &h_x, &h_y, &bed_elevation};
240 
241  unsigned int GHOSTS = eta.stencil_width();
242  assert(ice_thickness.stencil_width() >= GHOSTS);
243 
244  for (PointsWithGhosts p(*m_grid, GHOSTS); p; p.next()) {
245  const int i = p.i(), j = p.j();
246 
247  eta(i, j) = pow(ice_thickness(i, j), etapow);
248  }
249 
250  // now use Mahaffy on eta to get grad h on staggered;
251  // note grad h = (3/8) eta^{-5/8} grad eta + grad b because h = H + b
252 
253  assert(bed_elevation.stencil_width() >= 2);
254  assert(eta.stencil_width() >= 2);
255  assert(h_x.stencil_width() >= 1);
256  assert(h_y.stencil_width() >= 1);
257 
258  for (PointsWithGhosts p(*m_grid); p; p.next()) {
259  const int i = p.i(), j = p.j();
260 
261  auto b = bed_elevation.box(i, j);
262  auto e = eta.box(i, j);
263 
264  // i-offset
265  {
266  double mean_eta = 0.5 * (e.e + e.ij);
267  if (mean_eta > 0.0) {
268  double factor = invpow * pow(mean_eta, dinvpow);
269  h_x(i, j, 0) = factor * (e.e - e.ij) / dx;
270  h_y(i, j, 0) = factor * (e.ne + e.n - e.se - e.s) / (4.0 * dy);
271  } else {
272  h_x(i, j, 0) = 0.0;
273  h_y(i, j, 0) = 0.0;
274  }
275  // now add bed slope to get actual h_x, h_y
276  h_x(i, j, 0) += (b.e - b.ij) / dx;
277  h_y(i, j, 0) += (b.ne + b.n - b.se - b.s) / (4.0 * dy);
278  }
279 
280  // j-offset
281  {
282  double mean_eta = 0.5 * (e.n + e.ij);
283  if (mean_eta > 0.0) {
284  double factor = invpow * pow(mean_eta, dinvpow);
285  h_x(i, j, 1) = factor * (e.ne + e.e - e.nw - e.w) / (4.0 * dx);
286  h_y(i, j, 1) = factor * (e.n - e.ij) / dy;
287  } else {
288  h_x(i, j, 1) = 0.0;
289  h_y(i, j, 1) = 0.0;
290  }
291  // now add bed slope to get actual h_x, h_y
292  h_x(i, j, 1) += (b.ne + b.e - b.nw - b.w) / (4.0 * dx);
293  h_y(i, j, 1) += (b.n - b.ij) / dy;
294  }
295  } // end of the loop over grid points
296 }
297 
298 
299 //! \brief Compute the ice surface gradient using the Mary Anne Mahaffy method;
300 //! see [\ref Mahaffy].
301 void SIAFD::surface_gradient_mahaffy(const IceModelVec2S &ice_surface_elevation,
302  IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) {
303  const double dx = m_grid->dx(), dy = m_grid->dy(); // convenience
304 
305  const IceModelVec2S &h = ice_surface_elevation;
306 
307  IceModelVec::AccessList list{&h_x, &h_y, &h};
308 
309  // h_x and h_y have to have ghosts
310  assert(h_x.stencil_width() >= 1);
311  assert(h_y.stencil_width() >= 1);
312  // surface elevation needs more ghosts
313  assert(h.stencil_width() >= 2);
314 
315  for (PointsWithGhosts p(*m_grid); p; p.next()) {
316  const int i = p.i(), j = p.j();
317 
318  // I-offset
319  h_x(i, j, 0) = (h(i + 1, j) - h(i, j)) / dx;
320  h_y(i, j, 0) = (+ h(i + 1, j + 1) + h(i, j + 1)
321  - h(i + 1, j - 1) - h(i, j - 1)) / (4.0*dy);
322  // J-offset
323  h_y(i, j, 1) = (h(i, j + 1) - h(i, j)) / dy;
324  h_x(i, j, 1) = (+ h(i + 1, j + 1) + h(i + 1, j)
325  - h(i - 1, j + 1) - h(i - 1, j)) / (4.0*dx);
326  }
327 }
328 
329 //! \brief Compute the ice surface gradient using a modification of Marianne Haseloff's approach.
330 /*!
331  * The original code deals correctly with adjacent ice-free points with bed
332  * elevations which are above the surface of the ice nearby. This is done by
333  * setting surface gradient at the margin to zero at such locations.
334  *
335  * This code also deals with shelf fronts: sharp surface elevation change at
336  * the ice shelf front would otherwise cause abnormally high diffusivity
337  * values, which forces PISM to take shorter time-steps than necessary. (Note
338  * that the mass continuity code does not use SIA fluxes in floating areas.)
339  * This is done by assuming that the ice surface near shelf fronts is
340  * horizontal (i.e. here the surface gradient is set to zero also).
341  *
342  * The code below uses an interpretation of the standard Mahaffy scheme. We
343  * compute components of the surface gradient at staggered grid locations. The
344  * field h_x stores the x-component on the i-offset and j-offset grids, h_y ---
345  * the y-component.
346  *
347  * The Mahaffy scheme for the x-component at grid points on the i-offset grid
348  * (offset in the x-direction) is just the centered finite difference using
349  * adjacent regular-grid points. (Similarly for the y-component at j-offset
350  * locations.)
351  *
352  * Mahaffy's prescription for computing the y-component on the i-offset can be
353  * interpreted as:
354  *
355  * - compute the y-component at four surrounding j-offset staggered grid locations,
356  * - compute the average of these four.
357  *
358  * The code below does just that.
359  *
360  * - The first double for-loop computes x-components at i-offset
361  * locations and y-components at j-offset locations. Each computed
362  * number is assigned a weight (w_i and w_j) that is used below
363  *
364  * - The second double for-loop computes x-components at j-offset
365  * locations and y-components at i-offset locations as averages of
366  * quantities computed earlier. The weight are used to keep track of
367  * the number of values used in the averaging process.
368  *
369  * This method communicates ghost values of h_x and h_y. They cannot be
370  * computed locally because the first loop uses width=2 stencil of surface,
371  * mask, and bed to compute values at all grid points including width=1 ghosts,
372  * then the second loop uses width=1 stencil to compute local values. (In other
373  * words, a purely local computation would require width=3 stencil of surface,
374  * mask, and bed fields.)
375  */
376 void SIAFD::surface_gradient_haseloff(const IceModelVec2S &ice_surface_elevation,
377  const IceModelVec2CellType &cell_type,
378  IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) {
379  const double
380  dx = m_grid->dx(),
381  dy = m_grid->dy(); // convenience
382  const IceModelVec2S
383  &h = ice_surface_elevation;
385  &w_i = m_work_2d_0,
386  &w_j = m_work_2d_1; // averaging weights
387 
388  const IceModelVec2CellType &mask = cell_type;
389 
390  IceModelVec::AccessList list{&h_x, &h_y, &w_i, &w_j, &h, &mask};
391 
392  assert(mask.stencil_width() >= 2);
393  assert(h.stencil_width() >= 2);
394  assert(h_x.stencil_width() >= 1);
395  assert(h_y.stencil_width() >= 1);
396  assert(w_i.stencil_width() >= 1);
397  assert(w_j.stencil_width() >= 1);
398 
399  for (PointsWithGhosts p(*m_grid); p; p.next()) {
400  const int i = p.i(), j = p.j();
401 
402  // x-derivative, i-offset
403  {
404  if ((mask.floating_ice(i,j) && mask.ice_free_ocean(i+1,j)) ||
405  (mask.ice_free_ocean(i,j) && mask.floating_ice(i+1,j))) {
406  // marine margin
407  h_x(i,j,0) = 0;
408  w_i(i,j) = 0;
409  } else if ((mask.icy(i,j) && mask.ice_free(i+1,j) && h(i+1,j) > h(i,j)) ||
410  (mask.ice_free(i,j) && mask.icy(i+1,j) && h(i,j) > h(i+1,j))) {
411  // ice next to a "cliff"
412  h_x(i,j,0) = 0.0;
413  w_i(i,j) = 0;
414  } else {
415  // default case
416  h_x(i,j,0) = (h(i+1,j) - h(i,j)) / dx;
417  w_i(i,j) = 1;
418  }
419  }
420 
421  // y-derivative, j-offset
422  {
423  if ((mask.floating_ice(i,j) && mask.ice_free_ocean(i,j+1)) ||
424  (mask.ice_free_ocean(i,j) && mask.floating_ice(i,j+1))) {
425  // marine margin
426  h_y(i,j,1) = 0.0;
427  w_j(i,j) = 0.0;
428  } else if ((mask.icy(i,j) && mask.ice_free(i,j+1) && h(i,j+1) > h(i,j)) ||
429  (mask.ice_free(i,j) && mask.icy(i,j+1) && h(i,j) > h(i,j+1))) {
430  // ice next to a "cliff"
431  h_y(i,j,1) = 0.0;
432  w_j(i,j) = 0.0;
433  } else {
434  // default case
435  h_y(i,j,1) = (h(i,j+1) - h(i,j)) / dy;
436  w_j(i,j) = 1.0;
437  }
438  }
439  }
440 
441  for (Points p(*m_grid); p; p.next()) {
442  const int i = p.i(), j = p.j();
443 
444  // x-derivative, j-offset
445  {
446  if (w_j(i,j) > 0) {
447  double W = w_i(i,j) + w_i(i-1,j) + w_i(i-1,j+1) + w_i(i,j+1);
448  if (W > 0) {
449  h_x(i,j,1) = 1.0/W * (h_x(i,j,0) + h_x(i-1,j,0) + h_x(i-1,j+1,0) + h_x(i,j+1,0));
450  } else {
451  h_x(i,j,1) = 0.0;
452  }
453  } else {
454  if (mask.icy(i,j)) {
455  double W = w_i(i,j) + w_i(i-1,j);
456  if (W > 0) {
457  h_x(i,j,1) = 1.0/W * (h_x(i,j,0) + h_x(i-1,j,0));
458  } else {
459  h_x(i,j,1) = 0.0;
460  }
461  } else {
462  double W = w_i(i,j+1) + w_i(i-1,j+1);
463  if (W > 0) {
464  h_x(i,j,1) = 1.0/W * (h_x(i-1,j+1,0) + h_x(i,j+1,0));
465  } else {
466  h_x(i,j,1) = 0.0;
467  }
468  }
469  }
470  } // end of "x-derivative, j-offset"
471 
472  // y-derivative, i-offset
473  {
474  if (w_i(i,j) > 0) {
475  double W = w_j(i,j) + w_j(i,j-1) + w_j(i+1,j-1) + w_j(i+1,j);
476  if (W > 0) {
477  h_y(i,j,0) = 1.0/W * (h_y(i,j,1) + h_y(i,j-1,1) + h_y(i+1,j-1,1) + h_y(i+1,j,1));
478  } else {
479  h_y(i,j,0) = 0.0;
480  }
481  } else {
482  if (mask.icy(i,j)) {
483  double W = w_j(i,j) + w_j(i,j-1);
484  if (W > 0) {
485  h_y(i,j,0) = 1.0/W * (h_y(i,j,1) + h_y(i,j-1,1));
486  } else {
487  h_y(i,j,0) = 0.0;
488  }
489  } else {
490  double W = w_j(i+1,j-1) + w_j(i+1,j);
491  if (W > 0) {
492  h_y(i,j,0) = 1.0/W * (h_y(i+1,j-1,1) + h_y(i+1,j,1));
493  } else {
494  h_y(i,j,0) = 0.0;
495  }
496  }
497  }
498  } // end of "y-derivative, i-offset"
499  }
500 
501  h_x.update_ghosts();
502  h_y.update_ghosts();
503 }
504 
505 
506 //! \brief Compute the SIA diffusivity. If full_update, also store delta on the staggered grid.
507 /*!
508  * Recall that \f$ Q = -D \nabla h \f$ is the diffusive flux in the mass-continuity equation
509  *
510  * \f[ \frac{\partial H}{\partial t} = M - S - \nabla \cdot (Q + \mathbf{U}_b H),\f]
511  *
512  * where \f$h\f$ is the ice surface elevation, \f$M\f$ is the top surface
513  * accumulation/ablation rate, \f$S\f$ is the basal mass balance and
514  * \f$\mathbf{U}_b\f$ is the thickness-advective (in PISM: usually SSA) ice
515  * velocity.
516  *
517  * Recall also that at any particular point in the map-plane (i.e. if \f$x\f$
518  * and \f$y\f$ are fixed)
519  *
520  * \f[ D = 2\int_b^h F(z)P(z)(h-z)dz, \f]
521  *
522  * where \f$F(z)\f$ is a constitutive function and \f$P(z)\f$ is the pressure
523  * at a level \f$z\f$.
524  *
525  * By defining
526  *
527  * \f[ \delta(z) = 2F(z)P(z) \f]
528  *
529  * one can write
530  *
531  * \f[D = \int_b^h\delta(z)(h-z)dz. \f]
532  *
533  * The advantage is that it is then possible to avoid re-evaluating
534  * \f$F(z)\f$ (which is computationally expensive) in the horizontal ice
535  * velocity (see compute_3d_horizontal_velocity()) computation.
536  *
537  * This method computes \f$D\f$ and stores \f$\delta\f$ in delta[0,1] if full_update is true.
538  *
539  * The trapezoidal rule is used to approximate the integral.
540  *
541  * \param[in] full_update the flag specitying if we're doing a "full" update.
542  * \param[in] h_x x-component of the surface gradient, on the staggered grid
543  * \param[in] h_y y-component of the surface gradient, on the staggered grid
544  * \param[out] result diffusivity of the SIA flow
545  */
546 void SIAFD::compute_diffusivity(bool full_update,
547  const Geometry &geometry,
548  const IceModelVec3 *enthalpy,
549  const IceModelVec3 *age,
550  const IceModelVec2Stag &h_x,
551  const IceModelVec2Stag &h_y,
552  IceModelVec2Stag &result) {
554  &thk_smooth = m_work_2d_0,
555  &theta = m_work_2d_1;
556 
557  const IceModelVec2S
558  &h = geometry.ice_surface_elevation,
559  &H = geometry.ice_thickness;
560 
561  const IceModelVec2CellType &mask = geometry.cell_type;
562  IceModelVec3* delta[] = {&m_delta_0, &m_delta_1};
563 
564  result.set(0.0);
565 
566  const double
567  current_time = m_grid->ctx()->time()->current(),
568  D_limit = m_config->get_number("stress_balance.sia.max_diffusivity");
569 
570  const bool
571  compute_grain_size_using_age = m_config->get_flag("stress_balance.sia.grain_size_age_coupling"),
572  e_age_coupling = m_config->get_flag("stress_balance.sia.e_age_coupling"),
573  limit_diffusivity = m_config->get_flag("stress_balance.sia.limit_diffusivity"),
574  use_age = compute_grain_size_using_age or e_age_coupling;
575 
576  rheology::grain_size_vostok gs_vostok;
577 
578  // get "theta" from Schoof (2003) bed smoothness calculation and the
579  // thickness relative to the smoothed bed; each IceModelVec2S involved must
580  // have stencil width WIDE_GHOSTS for this too work
581  m_bed_smoother->theta(h, theta);
582 
583  m_bed_smoother->smoothed_thk(h, H, mask, thk_smooth);
584 
585  IceModelVec::AccessList list{&result, &theta, &thk_smooth, &h_x, &h_y, enthalpy};
586 
587  if (use_age) {
588  assert(age->stencil_width() >= 2);
589  list.add(*age);
590  }
591 
592  if (full_update) {
593  list.add({delta[0], delta[1]});
594  assert(m_delta_0.stencil_width() >= 1);
595  assert(m_delta_1.stencil_width() >= 1);
596  }
597 
598  assert(theta.stencil_width() >= 2);
599  assert(thk_smooth.stencil_width() >= 2);
600  assert(result.stencil_width() >= 1);
601  assert(h_x.stencil_width() >= 1);
602  assert(h_y.stencil_width() >= 1);
603  assert(enthalpy->stencil_width() >= 2);
604 
605  const std::vector<double> &z = m_grid->z();
606  const unsigned int
607  Mx = m_grid->Mx(),
608  My = m_grid->My(),
609  Mz = m_grid->Mz();
610 
611  std::vector<double> depth(Mz), stress(Mz), pressure(Mz), E(Mz), flow(Mz);
612  std::vector<double> delta_ij(Mz);
613  std::vector<double> A(Mz), ice_grain_size(Mz, m_config->get_number("constants.ice.grain_size", "m"));
614  std::vector<double> e_factor(Mz, m_e_factor);
615 
616  double D_max = 0.0;
617  int high_diffusivity_counter = 0;
618  for (int o=0; o<2; o++) {
619  ParallelSection loop(m_grid->com);
620  try {
621  for (PointsWithGhosts p(*m_grid); p; p.next()) {
622  const int i = p.i(), j = p.j();
623 
624  // staggered point: o=0 is i+1/2, o=1 is j+1/2, (i, j) and (i+oi, j+oj)
625  // are regular grid neighbors of a staggered point:
626  const int oi = 1 - o, oj = o;
627 
628  const double
629  thk = 0.5 * (thk_smooth(i, j) + thk_smooth(i+oi, j+oj));
630 
631  // zero thickness case:
632  if (thk == 0.0) {
633  result(i, j, o) = 0.0;
634  if (full_update) {
635  delta[o]->set_column(i, j, 0.0);
636  }
637  continue;
638  }
639 
640  const int ks = m_grid->kBelowHeight(thk);
641 
642  for (int k = 0; k <= ks; ++k) {
643  depth[k] = thk - z[k];
644  }
645 
646  // pressure added by the ice (i.e. pressure difference between the
647  // current level and the top of the column)
648  m_EC->pressure(depth, ks, pressure); // FIXME issue #15
649 
650  if (use_age) {
651  const double
652  *age_ij = age->get_column(i, j),
653  *age_offset = age->get_column(i+oi, j+oj);
654 
655  for (int k = 0; k <= ks; ++k) {
656  A[k] = 0.5 * (age_ij[k] + age_offset[k]);
657  }
658 
659  if (compute_grain_size_using_age) {
660  for (int k = 0; k <= ks; ++k) {
661  // convert age from seconds to years:
662  ice_grain_size[k] = gs_vostok(A[k] * m_seconds_per_year);
663  }
664  }
665 
666  if (e_age_coupling) {
667  for (int k = 0; k <= ks; ++k) {
668  const double accumulation_time = current_time - A[k];
669  if (interglacial(accumulation_time)) {
670  e_factor[k] = m_e_factor_interglacial;
671  } else {
672  e_factor[k] = m_e_factor;
673  }
674  }
675  }
676  }
677 
678  {
679  const double
680  *E_ij = enthalpy->get_column(i, j),
681  *E_offset = enthalpy->get_column(i+oi, j+oj);
682  for (int k = 0; k <= ks; ++k) {
683  E[k] = 0.5 * (E_ij[k] + E_offset[k]);
684  }
685  }
686 
687  const double alpha = sqrt(PetscSqr(h_x(i, j, o)) + PetscSqr(h_y(i, j, o)));
688  for (int k = 0; k <= ks; ++k) {
689  stress[k] = alpha * pressure[k];
690  }
691 
692  m_flow_law->flow_n(&stress[0], &E[0], &pressure[0], &ice_grain_size[0], ks + 1,
693  &flow[0]);
694 
695  const double theta_local = 0.5 * (theta(i, j) + theta(i+oi, j+oj));
696  for (int k = 0; k <= ks; ++k) {
697  delta_ij[k] = e_factor[k] * theta_local * 2.0 * pressure[k] * flow[k];
698  }
699 
700  double D = 0.0; // diffusivity for deformational SIA flow
701  {
702  for (int k = 1; k <= ks; ++k) {
703  // trapezoidal rule
704  const double dz = z[k] - z[k-1];
705  D += 0.5 * dz * ((depth[k] + dz) * delta_ij[k-1] + depth[k] * delta_ij[k]);
706  }
707  // finish off D with (1/2) dz (0 + (H-z[ks])*delta_ij[ks]), but dz=H-z[ks]:
708  const double dz = thk - z[ks];
709  D += 0.5 * dz * dz * delta_ij[ks];
710  }
711 
712  // Override diffusivity at the edges of the domain. (At these
713  // locations PISM uses ghost cells *beyond* the boundary of
714  // the computational domain. This does not matter if the ice
715  // does not extend all the way to the domain boundary, as in
716  // whole-ice-sheet simulations. In a regional setup, though,
717  // this adjustment lets us avoid taking very small time-steps
718  // because of the possible thickness and bed elevation
719  // "discontinuities" at the boundary.)
720  {
721  if ((i < 0 or i >= (int)Mx - 1) and
722  not (m_grid->periodicity() & X_PERIODIC)) {
723  D = 0.0;
724  }
725  if ((j < 0 or j >= (int)My - 1) and
726  not (m_grid->periodicity() & Y_PERIODIC)) {
727  D = 0.0;
728  }
729  }
730 
731  if (limit_diffusivity and D >= D_limit) {
732  D = D_limit;
733  high_diffusivity_counter += 1;
734  }
735 
736  D_max = std::max(D_max, D);
737 
738  result(i, j, o) = D;
739 
740  // if doing the full update, fill the delta column above the ice and
741  // store it:
742  if (full_update) {
743  for (unsigned int k = ks + 1; k < Mz; ++k) {
744  delta_ij[k] = 0.0;
745  }
746  delta[o]->set_column(i, j, &delta_ij[0]);
747  }
748  } // i, j-loop
749  } catch (...) {
750  loop.failed();
751  }
752  loop.check();
753  } // o-loop
754 
755  m_D_max = GlobalMax(m_grid->com, D_max);
756 
757  high_diffusivity_counter = GlobalSum(m_grid->com, high_diffusivity_counter);
758 
759  if (m_D_max > D_limit) {
760  // This can happen only if stress_balance.sia.limit_diffusivity is false (m_D_max <=
761  // D_limit when limiting is enabled).
762 
764  "Maximum diffusivity of SIA flow (%f m2/s) is too high.\n"
765  "This probably means that the bed elevation or the ice thickness is "
766  "too rough.\n"
767  "Increase stress_balance.sia.max_diffusivity to suppress this message.", m_D_max);
768 
769  }
770 
771  if (high_diffusivity_counter > 0) {
772  // This can happen only if stress_balance.sia.limit_diffusivity is true and this
773  // limiting mechanism was active (high_diffusivity_counter is incremented only if
774  // limit_diffusivity is true).
775 
776  m_log->message(2, " SIA diffusivity was capped at %.2f m2/s at %d locations.\n",
777  D_limit, high_diffusivity_counter);
778  }
779 }
780 
782  const IceModelVec2Stag &diffusivity,
783  IceModelVec2Stag &result) {
784 
785  IceModelVec::AccessList list{&diffusivity, &h_x, &h_y, &result};
786 
787  for (int o = 0; o < 2; o++) {
788  ParallelSection loop(m_grid->com);
789  try {
790  for (PointsWithGhosts p(*m_grid); p; p.next()) {
791  const int i = p.i(), j = p.j();
792 
793  const double slope = (o == 0) ? h_x(i, j, o) : h_y(i, j, o);
794 
795  result(i, j, o) = - diffusivity(i, j, o) * slope;
796  }
797  } catch (...) {
798  loop.failed();
799  }
800  loop.check();
801  } // o-loop
802 }
803 
804 //! \brief Compute I.
805 /*!
806  * This computes
807  * \f[ I(z) = \int_b^z\delta(s)ds.\f]
808  *
809  * Uses the trapezoidal rule to approximate the integral.
810  *
811  * See compute_diffusive_flux() for the definition of \f$\delta\f$.
812  *
813  * The result is stored in work_3d[0,1] and is used to compute the SIA component
814  * of the 3D-distributed horizontal ice velocity.
815  */
816 void SIAFD::compute_I(const Geometry &geometry) {
817 
818  IceModelVec2S &thk_smooth = m_work_2d_0;
820  IceModelVec3* delta[] = {&m_delta_0, &m_delta_1};
821 
822  const IceModelVec2S
823  &h = geometry.ice_surface_elevation,
824  &H = geometry.ice_thickness;
825 
826  const IceModelVec2CellType &mask = geometry.cell_type;
827 
828  m_bed_smoother->smoothed_thk(h, H, mask, thk_smooth);
829 
830  IceModelVec::AccessList list{delta[0], delta[1], I[0], I[1], &thk_smooth};
831 
832  assert(I[0]->stencil_width() >= 1);
833  assert(I[1]->stencil_width() >= 1);
834  assert(delta[0]->stencil_width() >= 1);
835  assert(delta[1]->stencil_width() >= 1);
836  assert(thk_smooth.stencil_width() >= 2);
837 
838  const unsigned int Mz = m_grid->Mz();
839 
840  std::vector<double> dz(Mz);
841  for (unsigned int k = 1; k < Mz; ++k) {
842  dz[k] = m_grid->z(k) - m_grid->z(k - 1);
843  }
844 
845  for (int o = 0; o < 2; ++o) {
846  ParallelSection loop(m_grid->com);
847  try {
848  for (PointsWithGhosts p(*m_grid); p; p.next()) {
849  const int i = p.i(), j = p.j();
850 
851  const int oi = 1 - o, oj = o;
852  const double
853  thk = 0.5 * (thk_smooth(i, j) + thk_smooth(i + oi, j + oj));
854 
855  const double *delta_ij = delta[o]->get_column(i, j);
856  double *I_ij = I[o]->get_column(i, j);
857 
858  const unsigned int ks = m_grid->kBelowHeight(thk);
859 
860  // within the ice:
861  I_ij[0] = 0.0;
862  double I_current = 0.0;
863  for (unsigned int k = 1; k <= ks; ++k) {
864  // trapezoidal rule
865  I_current += 0.5 * dz[k] * (delta_ij[k - 1] + delta_ij[k]);
866  I_ij[k] = I_current;
867  }
868 
869  // above the ice:
870  for (unsigned int k = ks + 1; k < Mz; ++k) {
871  I_ij[k] = I_current;
872  }
873  }
874  } catch (...) {
875  loop.failed();
876  }
877  loop.check();
878  } // o-loop
879 }
880 
881 //! \brief Compute horizontal components of the SIA velocity (in 3D).
882 /*!
883  * Recall that
884  *
885  * \f[ \mathbf{U}(z) = -2 \nabla h \int_b^z F(s)P(s)ds + \mathbf{U}_b,\f]
886  *
887  * which can be written in terms of \f$I(z)\f$ defined in compute_I():
888  *
889  * \f[ \mathbf{U}(z) = -I(z) \nabla h + \mathbf{U}_b. \f]
890  *
891  * \note This is one of the places where "hybridization" is done.
892  *
893  * \param[in] h_x the X-component of the surface gradient, on the staggered grid
894  * \param[in] h_y the Y-component of the surface gradient, on the staggered grid
895  * \param[in] sliding_velocity the thickness-advective velocity from the underlying stress balance module
896  * \param[out] u_out the X-component of the resulting horizontal velocity field
897  * \param[out] v_out the Y-component of the resulting horizontal velocity field
898  */
900  const IceModelVec2Stag &h_x,
901  const IceModelVec2Stag &h_y,
902  const IceModelVec2V &sliding_velocity,
903  IceModelVec3 &u_out, IceModelVec3 &v_out) {
904 
905  compute_I(geometry);
906  // after the compute_I() call work_3d[0,1] contains I on the staggered grid
908 
909  IceModelVec::AccessList list{&u_out, &v_out, &h_x, &h_y, &sliding_velocity, I[0], I[1]};
910 
911  const unsigned int Mz = m_grid->Mz();
912 
913  for (Points p(*m_grid); p; p.next()) {
914  const int i = p.i(), j = p.j();
915 
916  const double
917  *I_e = I[0]->get_column(i, j),
918  *I_w = I[0]->get_column(i - 1, j),
919  *I_n = I[1]->get_column(i, j),
920  *I_s = I[1]->get_column(i, j - 1);
921 
922  // Fetch values from 2D fields *outside* of the k-loop:
923  const double
924  h_x_w = h_x(i - 1, j, 0),
925  h_x_e = h_x(i, j, 0),
926  h_x_n = h_x(i, j, 1),
927  h_x_s = h_x(i, j - 1, 1);
928 
929  const double
930  h_y_w = h_y(i - 1, j, 0),
931  h_y_e = h_y(i, j, 0),
932  h_y_n = h_y(i, j, 1),
933  h_y_s = h_y(i, j - 1, 1);
934 
935  const double
936  sliding_velocity_u = sliding_velocity(i, j).u,
937  sliding_velocity_v = sliding_velocity(i, j).v;
938 
939  double
940  *u_ij = u_out.get_column(i, j),
941  *v_ij = v_out.get_column(i, j);
942 
943  // split into two loops to encourage auto-vectorization
944  for (unsigned int k = 0; k < Mz; ++k) {
945  u_ij[k] = sliding_velocity_u - 0.25 * (I_e[k] * h_x_e + I_w[k] * h_x_w +
946  I_n[k] * h_x_n + I_s[k] * h_x_s);
947  }
948  for (unsigned int k = 0; k < Mz; ++k) {
949  v_ij[k] = sliding_velocity_v - 0.25 * (I_e[k] * h_y_e + I_w[k] * h_y_w +
950  I_n[k] * h_y_n + I_s[k] * h_y_s);
951  }
952  }
953 
954  // Communicate to get ghosts:
955  u_out.update_ghosts();
956  v_out.update_ghosts();
957 }
958 
959 //! Determine if `accumulation_time` corresponds to an interglacial period.
960 bool SIAFD::interglacial(double accumulation_time) const {
961  if (accumulation_time < m_eemian_start) {
962  return false;
963  }
964 
965  if (accumulation_time < m_eemian_end) {
966  return true;
967  }
968 
969  return (accumulation_time >= m_holocene_start);
970 }
971 
973  return m_h_x;
974 }
975 
977  return m_h_y;
978 }
979 
981  return m_D;
982 }
983 
985  return *m_bed_smoother;
986 }
987 
988 
989 } // end of namespace stressbalance
990 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:140
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
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
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
bool floating_ice(int i, int j) const
bool ice_free(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
stencils::Box< double > box(int i, int j) const
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
double * get_column(int i, int j)
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition: iceModelVec3.cc:63
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:334
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void failed()
Indicates a failure of a parallel section.
void begin(const char *name) const
Definition: Profiling.cc:74
void end(const char *name) const
Definition: Profiling.cc:90
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
std::shared_ptr< FlowLaw > create()
A relationship between the age of the ice and the grain size from the Vostok core.
void preprocess_bed(const IceModelVec2S &topg)
Definition: BedSmoother.cc:99
void theta(const IceModelVec2S &usurf, IceModelVec2S &result) const
Definition: BedSmoother.cc:351
void smoothed_thk(const IceModelVec2S &usurf, const IceModelVec2S &thk, const IceModelVec2CellType &mask, IceModelVec2S &thksmooth) const
Computes a smoothed thickness map.
Definition: BedSmoother.cc:284
PISM bed smoother, plus bed roughness parameterization, based on Schoof (2003).
Definition: BedSmoother.hh:80
const IceModelVec3 * enthalpy
const IceModelVec3 * age
const IceModelVec2Stag & surface_gradient_x() const
Definition: SIAFD.cc:972
virtual void compute_surface_gradient(const Inputs &inputs, IceModelVec2Stag &h_x, IceModelVec2Stag &h_y)
Compute the ice surface gradient for the SIA.
Definition: SIAFD.cc:197
virtual void surface_gradient_mahaffy(const IceModelVec2S &ice_surface_elevation, IceModelVec2Stag &h_x, IceModelVec2Stag &h_y)
Compute the ice surface gradient using the Mary Anne Mahaffy method; see [Mahaffy].
Definition: SIAFD.cc:301
virtual void compute_3d_horizontal_velocity(const Geometry &geometry, const IceModelVec2Stag &h_x, const IceModelVec2Stag &h_y, const IceModelVec2V &vel_input, IceModelVec3 &u_out, IceModelVec3 &v_out)
Compute horizontal components of the SIA velocity (in 3D).
Definition: SIAFD.cc:899
const IceModelVec2Stag & diffusivity() const
Definition: SIAFD.cc:980
const unsigned int m_stencil_width
Definition: SIAFD.hh:106
virtual void update(const IceModelVec2V &sliding_velocity, const Inputs &inputs, bool full_update)
Do the update; if full_update == false skip the update of 3D velocities and strain heating.
Definition: SIAFD.cc:125
bool interglacial(double accumulation_time) const
Determine if accumulation_time corresponds to an interglacial period.
Definition: SIAFD.cc:960
IceModelVec2Stag m_h_y
Definition: SIAFD.hh:112
virtual void surface_gradient_eta(const IceModelVec2S &ice_thickness, const IceModelVec2S &bed_elevation, IceModelVec2Stag &h_x, IceModelVec2Stag &h_y)
Compute the ice surface gradient using the eta-transformation.
Definition: SIAFD.cc:227
IceModelVec2S m_work_2d_1
Definition: SIAFD.hh:110
const BedSmoother & bed_smoother() const
Definition: SIAFD.cc:984
IceModelVec2Stag m_D
Definition: SIAFD.hh:112
IceModelVec2Stag m_h_x
temporary storage for the surface gradient and the diffusivity
Definition: SIAFD.hh:112
virtual void surface_gradient_haseloff(const IceModelVec2S &ice_surface_elevation, const IceModelVec2CellType &cell_type, IceModelVec2Stag &h_x, IceModelVec2Stag &h_y)
Compute the ice surface gradient using a modification of Marianne Haseloff's approach.
Definition: SIAFD.cc:376
IceModelVec3 m_delta_1
Definition: SIAFD.hh:115
double m_e_factor_interglacial
Definition: SIAFD.hh:133
BedSmoother * m_bed_smoother
Definition: SIAFD.hh:120
IceModelVec2S m_work_2d_0
temporary storage for eta, theta and the smoothed thickness
Definition: SIAFD.hh:109
virtual void init()
Initialize the SIA module.
Definition: SIAFD.cc:102
virtual void compute_diffusivity(bool full_update, const Geometry &geometry, const IceModelVec3 *enthalpy, const IceModelVec3 *age, const IceModelVec2Stag &h_x, const IceModelVec2Stag &h_y, IceModelVec2Stag &result)
Compute the SIA diffusivity. If full_update, also store delta on the staggered grid.
Definition: SIAFD.cc:546
SIAFD(IceGrid::ConstPtr g)
Definition: SIAFD.cc:43
IceModelVec3 m_work_3d_0
temporary storage used to store I and strain_heating on the staggered grid
Definition: SIAFD.hh:117
virtual void compute_I(const Geometry &geometry)
Compute I.
Definition: SIAFD.cc:816
IceModelVec3 m_delta_0
temporary storage for delta on the staggered grid
Definition: SIAFD.hh:114
virtual void compute_diffusive_flux(const IceModelVec2Stag &h_x, const IceModelVec2Stag &h_y, const IceModelVec2Stag &diffusivity, IceModelVec2Stag &result)
Definition: SIAFD.cc:781
IceModelVec3 m_work_3d_1
Definition: SIAFD.hh:118
const IceModelVec2Stag & surface_gradient_y() const
Definition: SIAFD.cc:976
EnthalpyConverter::Ptr m_EC
Definition: SSB_Modifier.hh:66
std::shared_ptr< rheology::FlowLaw > m_flow_law
Definition: SSB_Modifier.hh:65
Shallow stress balance modifier (such as the non-sliding SIA).
Definition: SSB_Modifier.hh:39
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
bool FlowLawUsesGrainSize(const FlowLaw &flow_law)
Definition: FlowLaw.cc:263
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:72
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
@ Y_PERIODIC
Definition: IceGrid.hh:48
@ X_PERIODIC
Definition: IceGrid.hh:48
static const double g
Definition: exactTestP.cc:39
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
static const double k
Definition: exactTestP.cc:45
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
const int I[]
Definition: ssafd_code.cc:24