PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
StressBalance.cc
Go to the documentation of this file.
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022 Constantine Khroulev and Ed Bueler
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 "StressBalance.hh"
20 #include "ShallowStressBalance.hh"
21 #include "SSB_Modifier.hh"
22 #include "pism/coupler/OceanModel.hh"
23 #include "pism/util/EnthalpyConverter.hh"
24 #include "pism/rheology/FlowLaw.hh"
25 #include "pism/util/IceGrid.hh"
26 #include "pism/util/Mask.hh"
27 #include "pism/util/ConfigInterface.hh"
28 #include "pism/util/Vars.hh"
29 #include "pism/util/error_handling.hh"
30 #include "pism/util/Profiling.hh"
31 #include "pism/util/IceModelVec2CellType.hh"
32 #include "pism/util/Time.hh"
33 #include "pism/geometry/Geometry.hh"
34 #include "pism/util/Context.hh"
35 
36 namespace pism {
37 namespace stressbalance {
38 
40  geometry = NULL;
41  new_bed_elevation = true;
42 
43  basal_melt_rate = NULL;
44  water_column_pressure = NULL;
45  fracture_density = NULL;
46  basal_yield_stress = NULL;
47 
48  enthalpy = NULL;
49  age = NULL;
50 
51  bc_mask = NULL;
52  bc_values = NULL;
53 
54  no_model_mask = NULL;
57 }
58 
59 /*!
60  * Save stress balance inputs to a file (for debugging).
61  */
62 void Inputs::dump(const char *filename) const {
63  if (not geometry) {
64  return;
65  }
66 
67  auto ctx = geometry->ice_thickness.grid()->ctx();
68  auto config = ctx->config();
69 
70  File output(ctx->com(), filename,
71  string_to_backend(config->get_string("output.format")),
73 
74  config->write(output);
75 
76  io::define_time(output, *ctx);
77  io::append_time(output, config->get_string("time.dimension_name"), ctx->time()->current());
78 
79  {
80  geometry->latitude.write(output);
81  geometry->longitude.write(output);
82 
83  geometry->bed_elevation.write(output);
85 
86  geometry->ice_thickness.write(output);
88 
89  geometry->cell_type.write(output);
92  }
93 
94  if (basal_melt_rate) {
95  basal_melt_rate->write(output);
96  }
97 
100  }
101 
102  if (fracture_density) {
103  fracture_density->write(output);
104  }
105 
106  if (basal_yield_stress) {
107  basal_yield_stress->write(output);
108  }
109 
110  if (enthalpy) {
111  enthalpy->write(output);
112  }
113 
114  if (age) {
115  age->write(output);
116  }
117 
118  if (bc_mask) {
119  bc_mask->write(output);
120  }
121 
122  if (bc_values) {
123  bc_values->write(output);
124  }
125 
126  if (no_model_mask) {
127  no_model_mask->write(output);
128  }
129 
131  no_model_ice_thickness->write(output);
132  }
133 
136  }
137 }
138 
140  std::shared_ptr<ShallowStressBalance> sb,
141  std::shared_ptr<SSB_Modifier> ssb_mod)
142  : Component(g),
143  m_w(m_grid, "wvel_rel", WITHOUT_GHOSTS, m_grid->z()),
144  m_strain_heating(m_grid, "strain_heating", WITHOUT_GHOSTS, m_grid->z()),
145  m_shallow_stress_balance(sb),
146  m_modifier(ssb_mod) {
147 
148  m_w.set_attrs("diagnostic",
149  "vertical velocity of ice, relative to base of ice directly below",
150  "m s-1", "m year-1", "", 0);
151  m_w.set_time_independent(false);
152 
153  m_strain_heating.set_attrs("internal",
154  "rate of strain heating in ice (dissipation heating)",
155  "W m-3", "W m-3", "", 0);
156 }
157 
159 }
160 
161 //! \brief Initialize the StressBalance object.
163  m_shallow_stress_balance->init();
164  m_modifier->init();
165 }
166 
167 //! \brief Performs the shallow stress balance computation.
168 void StressBalance::update(const Inputs &inputs, bool full_update) {
169 
170  const Profiling &profiling = m_grid->ctx()->profiling();
171 
172  try {
173  profiling.begin("stress_balance.shallow");
174  m_shallow_stress_balance->update(inputs, full_update);
175  profiling.end("stress_balance.shallow");
176 
177  profiling.begin("stress_balance.modifier");
178  m_modifier->update(m_shallow_stress_balance->velocity(),
179  inputs, full_update);
180  profiling.end("stress_balance.modifier");
181 
182  if (full_update) {
183  const IceModelVec3 &u = m_modifier->velocity_u();
184  const IceModelVec3 &v = m_modifier->velocity_v();
185 
186  profiling.begin("stress_balance.strain_heat");
187  this->compute_volumetric_strain_heating(inputs);
188  profiling.end("stress_balance.strain_heat");
189 
190  profiling.begin("stress_balance.vertical_velocity");
192  u, v, inputs.basal_melt_rate, m_w);
193  profiling.end("stress_balance.vertical_velocity");
194 
196  inputs.geometry->cell_type,
197  u, v, m_w);
198  }
199 
201  inputs.geometry->cell_type,
202  m_shallow_stress_balance->velocity());
203  }
204  catch (RuntimeError &e) {
205  e.add_context("updating the stress balance");
206  throw;
207  }
208 }
209 
211  return m_cfl_2d;
212 }
213 
215  return m_cfl_3d;
216 }
217 
219  return m_shallow_stress_balance->velocity();
220 }
221 
223  return m_modifier->diffusive_flux();
224 }
225 
227  return m_modifier->max_diffusivity();
228 }
229 
231  return m_modifier->velocity_u();
232 }
233 
235  return m_modifier->velocity_v();
236 }
237 
239  return m_w;
240 }
241 
243  return m_shallow_stress_balance->basal_frictional_heating();
244 }
245 
247  return m_strain_heating;
248 }
249 
250 //! Compute vertical velocity using incompressibility of the ice.
251 /*!
252 The vertical velocity \f$w(x,y,z,t)\f$ is the velocity *relative to the
253 location of the base of the ice column*. That is, the vertical velocity
254 computed here is identified as \f$\tilde w(x,y,s,t)\f$ in the page
255 []@ref vertchange.
256 
257 Thus \f$w<0\f$ here means that that
258 that part of the ice is getting closer to the base of the ice, and so on.
259 The slope of the bed (i.e. relative to the geoid) and/or the motion of the
260 bed (i.e. from bed deformation) do not affect the vertical velocity.
261 
262 In fact the following statement is exactly true if the basal melt rate is zero:
263 the vertical velocity at a point in the ice is positive (negative) if and only
264 if the average horizontal divergence of the horizontal velocity, in the portion
265 of the ice column below that point, is negative (positive).
266 In particular, because \f$z=0\f$ is the location of the base of the ice
267 always, the only way to have \f$w(x,y,0,t) \ne 0\f$ is to have a basal melt
268 rate.
269 
270 Incompressibility itself says
271  \f[ \nabla\cdot\mathbf{U} + \frac{\partial w}{\partial z} = 0. \f]
272 This is immediately equivalent to the integral
273  \f[ w(x,y,z,t) = - \int_{b(x,y,t)}^{z} \nabla\cdot\mathbf{U}\,d\zeta
274  + w_b(x,y,t). \f]
275 Here the value \f$w_b(x,y,t)\f$ is either zero or the negative of the basal melt rate
276 according to the value of the flag `geometry.update.use_basal_melt_rate`.
277 
278 The vertical integral is computed by the trapezoid rule.
279  */
281  const IceModelVec3 &u,
282  const IceModelVec3 &v,
283  const IceModelVec2S *basal_melt_rate,
284  IceModelVec3 &result) {
285 
286  const bool use_upstream_fd = m_config->get_string("stress_balance.vertical_velocity_approximation") == "upstream";
287 
288  IceModelVec::AccessList list{&u, &v, &mask, &result};
289 
290  if (basal_melt_rate) {
291  list.add(*basal_melt_rate);
292  }
293 
294  const std::vector<double> &z = m_grid->z();
295  const unsigned int Mz = m_grid->Mz();
296 
297  const double
298  dx = m_grid->dx(),
299  dy = m_grid->dy();
300 
301  std::vector<double> u_x_plus_v_y(Mz);
302 
303  for (Points p(*m_grid); p; p.next()) {
304  const int i = p.i(), j = p.j();
305 
306  double *w_ij = result.get_column(i,j);
307 
308  const double
309  *u_w = u.get_column(i-1,j),
310  *u_ij = u.get_column(i,j),
311  *u_e = u.get_column(i+1,j);
312  const double
313  *v_s = v.get_column(i,j-1),
314  *v_ij = v.get_column(i,j),
315  *v_n = v.get_column(i,j+1);
316 
317  double
318  west = 1.0,
319  east = 1.0,
320  south = 1.0,
321  north = 1.0;
322  double
323  D_x = 0, // 1/(dx), 1/(2dx), or 0
324  D_y = 0; // 1/(dy), 1/(2dy), or 0
325 
326  // Switch between second-order centered differences in the interior and
327  // first-order one-sided differences at ice margins.
328 
329  // x-derivative
330  {
331  // use basal velocity to determine FD direction ("upwind" when it's clear, centered when it's
332  // not)
333  if (use_upstream_fd) {
334  const double
335  uw = 0.5 * (u_w[0] + u_ij[0]),
336  ue = 0.5 * (u_ij[0] + u_e[0]);
337 
338  if (uw > 0.0 and ue >= 0.0) {
339  west = 1.0;
340  east = 0.0;
341  } else if (uw <= 0.0 and ue < 0.0) {
342  west = 0.0;
343  east = 1.0;
344  } else {
345  west = 1.0;
346  east = 1.0;
347  }
348  }
349 
350  if ((mask.icy(i,j) and mask.ice_free(i+1,j)) or (mask.ice_free(i,j) and mask.icy(i+1,j))) {
351  east = 0;
352  }
353  if ((mask.icy(i,j) and mask.ice_free(i-1,j)) or (mask.ice_free(i,j) and mask.icy(i-1,j))) {
354  west = 0;
355  }
356 
357  if (east + west > 0) {
358  D_x = 1.0 / (dx * (east + west));
359  } else {
360  D_x = 0.0;
361  }
362  }
363 
364  // y-derivative
365  {
366  // use basal velocity to determine FD direction ("upwind" when it's clear, centered when it's
367  // not)
368  if (use_upstream_fd) {
369  const double
370  vs = 0.5 * (v_s[0] + v_ij[0]),
371  vn = 0.5 * (v_ij[0] + v_n[0]);
372 
373  if (vs > 0.0 and vn >= 0.0) {
374  south = 1.0;
375  north = 0.0;
376  } else if (vs <= 0.0 and vn < 0.0) {
377  south = 0.0;
378  north = 1.0;
379  } else {
380  south = 1.0;
381  north = 1.0;
382  }
383  }
384 
385  if ((mask.icy(i,j) and mask.ice_free(i,j+1)) or (mask.ice_free(i,j) and mask.icy(i,j+1))) {
386  north = 0;
387  }
388  if ((mask.icy(i,j) and mask.ice_free(i,j-1)) or (mask.ice_free(i,j) and mask.icy(i,j-1))) {
389  south = 0;
390  }
391 
392  if (north + south > 0) {
393  D_y = 1.0 / (dy * (north + south));
394  } else {
395  D_y = 0.0;
396  }
397  }
398 
399  // compute u_x + v_y using a vectorizable loop
400  for (unsigned int k = 0; k < Mz; ++k) {
401  double
402  u_x = D_x * (west * (u_ij[k] - u_w[k]) + east * (u_e[k] - u_ij[k])),
403  v_y = D_y * (south * (v_ij[k] - v_s[k]) + north * (v_n[k] - v_ij[k]));
404  u_x_plus_v_y[k] = u_x + v_y;
405  }
406 
407  // at the base: include the basal melt rate
408  if (basal_melt_rate != NULL) {
409  w_ij[0] = - (*basal_melt_rate)(i,j);
410  } else {
411  w_ij[0] = 0.0;
412  }
413 
414  // within the ice and above:
415  for (unsigned int k = 1; k < Mz; ++k) {
416  const double dz = z[k] - z[k-1];
417 
418  w_ij[k] = w_ij[k - 1] - (0.5 * dz) * (u_x_plus_v_y[k] + u_x_plus_v_y[k - 1]);
419  }
420  }
421 }
422 
423 /**
424  * This function computes \f$D^2\f$ defined by
425  *
426  * \f[ 2D^2 = D_{ij} D_{ij}\f]
427  * or
428  * \f[
429  * D^2 = \frac{1}{2}\,\left(\frac{1}{2}\,(v_{z})^2 + (v_{y} + u_{x})^2 +
430  * (v_{y})^2 + \frac{1}{2}\,(v_{x} + u_{y})^2 + \frac{1}{2}\,(u_{z})^2 +
431  * (u_{x})^2\right)
432  * \f]
433  *
434  * (note the use of the summation convention). Here \f$D_{ij}\f$ is the
435  * strain rate tensor. See
436  * StressBalance::compute_volumetric_strain_heating() for details.
437  *
438  * @param u_x,u_y,u_z partial derivatives of \f$u\f$, the x-component of the ice velocity
439  * @param v_x,v_y,v_z partial derivatives of \f$v\f$, the y-component of the ice velocity
440  *
441  * @return \f$D^2\f$, where \f$D\f$ is defined above.
442  */
443 static inline double D2(double u_x, double u_y, double u_z, double v_x, double v_y, double v_z) {
444  return 0.5 * (PetscSqr(u_x + v_y) + u_x*u_x + v_y*v_y + 0.5 * (PetscSqr(u_y + v_x) + u_z*u_z + v_z*v_z));
445 }
446 
447 /**
448  \brief Computes the volumetric strain heating using horizontal
449  velocity.
450 
451  Following the notation used in [\ref BBssasliding], let \f$u\f$ be a
452  three-dimensional *vector* velocity field. Then the strain rate
453  tensor \f$D_{ij}\f$ is defined by
454 
455  \f[ D_{ij} = \frac 12 \left(\diff{u_{i}}{x_{j}} + \diff{u_{j}}{x_{i}} \right), \f]
456 
457  Where \f$i\f$ and \f$j\f$ range from \f$1\f$ to \f$3\f$.
458 
459  The flow law in the viscosity form states
460 
461  \f[ \tau_{ij} = 2 \eta D_{ij}, \f]
462 
463  and the nonlinear ice viscosity satisfies
464 
465  \f[ 2 \eta = B(T) D^{(1/n) - 1}. \f]
466 
467  Here \f$D^{2}\f$ is defined by \f$2D^{2} = D_{ij}D_{ij}\f$ (using the
468  summation convention) and \f$B(T) = A(T)^{-1/n}\f$ is the ice hardness.
469 
470  Now the volumetric strain heating is
471 
472  \f[ \Sigma = \sum_{i,j=1}^{3}D_{ij}\tau_{ij} = 2 B(T) D^{(1/n) + 1}. \f]
473 
474  We use an *approximation* of \f$D_{ij}\f$ common in shallow ice models:
475 
476  - we assume that horizontal derivatives of the vertical velocity are
477  much smaller than \f$z\f$ derivatives horizontal velocity
478  components \f$u\f$ and \f$v\f$. (We drop \f$w_x\f$ and \f$w_y\f$
479  terms in \f$D_{ij}\f$.)
480 
481  - we use the incompressibility of ice to approximate \f$w_z\f$:
482 
483  \f[ w_z = - (u_x + v_y). \f]
484 
485  Requires ghosts of `u` and `v` velocity components and uses the fact
486  that `u` and `v` above the ice are filled using constant
487  extrapolation.
488 
489  Resulting field does not have ghosts.
490 
491  Below is the *Maxima* code that produces the expression evaluated by D2().
492 
493  derivabbrev : true;
494  U : [u, v, w]; X : [x, y, z]; depends(U, X);
495  gradef(w, x, 0); gradef(w, y, 0);
496  gradef(w, z, -(diff(u, x) + diff(v, y)));
497  d[i,j] := 1/2 * (diff(U[i], X[j]) + diff(U[j], X[i]));
498  D : genmatrix(d, 3, 3), ratsimp, factor;
499  tex('D = D);
500  tex('D^2 = 1/2 * mat_trace(D . D));
501 
502  @return 0 on success
503  */
505  PetscErrorCode ierr;
506 
507  const rheology::FlowLaw &flow_law = *m_shallow_stress_balance->flow_law();
508  EnthalpyConverter::Ptr EC = m_shallow_stress_balance->enthalpy_converter();
509 
510  const IceModelVec3
511  &u = m_modifier->velocity_u(),
512  &v = m_modifier->velocity_v();
513 
514  const IceModelVec2S &thickness = inputs.geometry->ice_thickness;
515  const IceModelVec3 *enthalpy = inputs.enthalpy;
516 
517  const IceModelVec2CellType &mask = inputs.geometry->cell_type;
518 
519  double
520  enhancement_factor = m_shallow_stress_balance->flow_enhancement_factor(),
521  n = flow_law.exponent(),
522  exponent = 0.5 * (1.0 / n + 1.0),
523  e_to_a_power = pow(enhancement_factor,-1.0/n);
524 
525  IceModelVec::AccessList list{&mask, enthalpy, &m_strain_heating, &thickness, &u, &v};
526 
527  const std::vector<double> &z = m_grid->z();
528  const unsigned int Mz = m_grid->Mz();
529  std::vector<double> depth(Mz), pressure(Mz), hardness(Mz);
530 
531  ParallelSection loop(m_grid->com);
532  try {
533  for (Points p(*m_grid); p; p.next()) {
534  const int i = p.i(), j = p.j();
535 
536  double H = thickness(i, j);
537  int ks = m_grid->kBelowHeight(H);
538  const double
539  *u_ij, *u_w, *u_n, *u_e, *u_s,
540  *v_ij, *v_w, *v_n, *v_e, *v_s;
541  double *Sigma;
542  const double *E_ij;
543 
544  double west = 1, east = 1, south = 1, north = 1,
545  D_x = 0, // 1/(dx), 1/(2dx), or 0
546  D_y = 0; // 1/(dy), 1/(2dy), or 0
547 
548  // x-derivative
549  {
550  if ((mask.icy(i,j) and mask.ice_free(i+1,j)) or (mask.ice_free(i,j) and mask.icy(i+1,j))) {
551  east = 0;
552  }
553  if ((mask.icy(i,j) and mask.ice_free(i-1,j)) or (mask.ice_free(i,j) and mask.icy(i-1,j))) {
554  west = 0;
555  }
556 
557  if (east + west > 0) {
558  D_x = 1.0 / (m_grid->dx() * (east + west));
559  } else {
560  D_x = 0.0;
561  }
562  }
563 
564  // y-derivative
565  {
566  if ((mask.icy(i,j) and mask.ice_free(i,j+1)) or (mask.ice_free(i,j) and mask.icy(i,j+1))) {
567  north = 0;
568  }
569  if ((mask.icy(i,j) and mask.ice_free(i,j-1)) or (mask.ice_free(i,j) and mask.icy(i,j-1))) {
570  south = 0;
571  }
572 
573  if (north + south > 0) {
574  D_y = 1.0 / (m_grid->dy() * (north + south));
575  } else {
576  D_y = 0.0;
577  }
578  }
579 
580  u_ij = u.get_column(i, j);
581  u_w = u.get_column(i - 1, j);
582  u_e = u.get_column(i + 1, j);
583  u_s = u.get_column(i, j - 1);
584  u_n = u.get_column(i, j + 1);
585 
586  v_ij = v.get_column(i, j);
587  v_w = v.get_column(i - 1, j);
588  v_e = v.get_column(i + 1, j);
589  v_s = v.get_column(i, j - 1);
590  v_n = v.get_column(i, j + 1);
591 
592  E_ij = enthalpy->get_column(i, j);
593  Sigma = m_strain_heating.get_column(i, j);
594 
595  for (int k = 0; k <= ks; ++k) {
596  depth[k] = H - z[k];
597  }
598 
599  // pressure added by the ice (i.e. pressure difference between the
600  // current level and the top of the column)
601  EC->pressure(depth, ks, pressure); // FIXME issue #15
602 
603  flow_law.hardness_n(E_ij, &pressure[0], ks + 1, &hardness[0]);
604 
605  for (int k = 0; k <= ks; ++k) {
606  double dz;
607 
608  double u_z = 0.0, v_z = 0.0,
609  u_x = D_x * (west * (u_ij[k] - u_w[k]) + east * (u_e[k] - u_ij[k])),
610  u_y = D_y * (south * (u_ij[k] - u_s[k]) + north * (u_n[k] - u_ij[k])),
611  v_x = D_x * (west * (v_ij[k] - v_w[k]) + east * (v_e[k] - v_ij[k])),
612  v_y = D_y * (south * (v_ij[k] - v_s[k]) + north * (v_n[k] - v_ij[k]));
613 
614  if (k > 0) {
615  dz = z[k+1] - z[k-1];
616  u_z = (u_ij[k+1] - u_ij[k-1]) / dz;
617  v_z = (v_ij[k+1] - v_ij[k-1]) / dz;
618  } else {
619  // use one-sided differences for u_z and v_z on the bottom level
620  dz = z[1] - z[0];
621  u_z = (u_ij[1] - u_ij[0]) / dz;
622  v_z = (v_ij[1] - v_ij[0]) / dz;
623  }
624 
625  Sigma[k] = 2.0 * e_to_a_power * hardness[k] * pow(D2(u_x, u_y, u_z, v_x, v_y, v_z), exponent);
626  } // k-loop
627 
628  int remaining_levels = Mz - (ks + 1);
629  if (remaining_levels > 0) {
630 #if PETSC_VERSION_LT(3, 12, 0)
631  ierr = PetscMemzero(&Sigma[ks+1],
632  remaining_levels*sizeof(double));
633 #else
634  ierr = PetscArrayzero(&Sigma[ks+1], remaining_levels);
635 #endif
636  PISM_CHK(ierr, "PetscMemzero");
637  }
638  }
639  } catch (...) {
640  loop.failed();
641  }
642  loop.check();
643 }
644 
645 std::string StressBalance::stdout_report() const {
646  return m_shallow_stress_balance->stdout_report() + m_modifier->stdout_report();
647 }
648 
650  return m_shallow_stress_balance.get();
651 }
652 
654  return m_modifier.get();
655 }
656 
657 
658 void StressBalance::define_model_state_impl(const File &output) const {
659  m_shallow_stress_balance->define_model_state(output);
660  m_modifier->define_model_state(output);
661 }
662 
663 void StressBalance::write_model_state_impl(const File &output) const {
664  m_shallow_stress_balance->write_model_state(output);
665  m_modifier->write_model_state(output);
666 }
667 
668 //! \brief Compute eigenvalues of the horizontal, vertically-integrated strain rate tensor.
669 /*!
670 Calculates all components \f$D_{xx}, D_{yy}, D_{xy}=D_{yx}\f$ of the
671 vertically-averaged strain rate tensor \f$D\f$ [\ref SchoofStream]. Then computes
672 the eigenvalues `result(i,j,0)` = (maximum eigenvalue), `result(i,j,1)` = (minimum
673 eigenvalue). Uses the provided thickness to make decisions (PIK) about computing
674 strain rates near calving front.
675 
676 Note that `result(i,j,0)` >= `result(i,j,1)`, but there is no necessary relation between
677 the magnitudes, and either principal strain rate could be negative or positive.
678 
679 Result can be used in a calving law, for example in eigencalving (PIK).
680 
681 Note: strain rates will be derived from SSA velocities, using ghosts when
682 necessary. Both implementations (SSAFD and SSAFEM) call
683 update_ghosts() to ensure that ghost values are up to date.
684  */
686  const IceModelVec2CellType &mask,
687  IceModelVec3 &result) {
688 
689  using mask::ice_free;
690 
691  IceGrid::ConstPtr grid = result.grid();
692  double dx = grid->dx();
693  double dy = grid->dy();
694 
695  if (result.ndof() != 2) {
696  throw RuntimeError(PISM_ERROR_LOCATION, "result.dof() == 2 is required");
697  }
698 
699  IceModelVec::AccessList list{&V, &mask, &result};
700 
701  for (Points p(*grid); p; p.next()) {
702  const int i = p.i(), j = p.j();
703 
704  if (mask.ice_free(i,j)) {
705  result(i,j,0) = 0.0;
706  result(i,j,1) = 0.0;
707  continue;
708  }
709 
710  auto m = mask.star(i,j);
711  auto U = V.star(i,j);
712 
713  // strain in units s-1
714  double u_x = 0, u_y = 0, v_x = 0, v_y = 0,
715  east = 1, west = 1, south = 1, north = 1;
716 
717  // Computes u_x using second-order centered finite differences written as
718  // weighted sums of first-order one-sided finite differences.
719  //
720  // Given the cell layout
721  // *----n----*
722  // | |
723  // | |
724  // w e
725  // | |
726  // | |
727  // *----s----*
728  // east == 0 if the east neighbor of the current cell is ice-free. In
729  // this case we use the left- (west-) sided difference.
730  //
731  // If both neighbors in the east-west (x) direction are ice-free the
732  // x-derivative is set to zero (see u_x, v_x initialization above).
733  //
734  // Similarly in other directions.
735  if (ice_free(m.e)) {
736  east = 0;
737  }
738  if (ice_free(m.w)) {
739  west = 0;
740  }
741  if (ice_free(m.n)) {
742  north = 0;
743  }
744  if (ice_free(m.s)) {
745  south = 0;
746  }
747 
748  if (west + east > 0) {
749  u_x = 1.0 / (dx * (west + east)) * (west * (U.ij.u - U[West].u) + east * (U[East].u - U.ij.u));
750  v_x = 1.0 / (dx * (west + east)) * (west * (U.ij.v - U[West].v) + east * (U[East].v - U.ij.v));
751  }
752 
753  if (south + north > 0) {
754  u_y = 1.0 / (dy * (south + north)) * (south * (U.ij.u - U[South].u) + north * (U[North].u - U.ij.u));
755  v_y = 1.0 / (dy * (south + north)) * (south * (U.ij.v - U[South].v) + north * (U[North].v - U.ij.v));
756  }
757 
758  const double A = 0.5 * (u_x + v_y), // A = (1/2) trace(D)
759  B = 0.5 * (u_x - v_y),
760  Dxy = 0.5 * (v_x + u_y), // B^2 = A^2 - u_x v_y
761  q = sqrt(B*B + Dxy*Dxy);
762  result(i,j,0) = A + q;
763  result(i,j,1) = A - q; // q >= 0 so e1 >= e2
764 
765  }
766 }
767 
768 //! @brief Compute 2D deviatoric stresses.
769 /*! Note: IceModelVec2 result has to have dof == 3. */
771  const IceModelVec2V &velocity,
772  const IceModelVec2S &hardness,
773  const IceModelVec2CellType &cell_type,
774  IceModelVec3 &result) {
775 
776  using mask::ice_free;
777 
778  auto grid = result.grid();
779 
780  const double
781  dx = grid->dx(),
782  dy = grid->dy();
783 
784  if (result.ndof() != 3) {
785  throw RuntimeError(PISM_ERROR_LOCATION, "result.ndof() == 3 is required");
786  }
787 
788  IceModelVec::AccessList list{&velocity, &hardness, &result, &cell_type};
789 
790  for (Points p(*grid); p; p.next()) {
791  const int i = p.i(), j = p.j();
792 
793  if (cell_type.ice_free(i, j)) {
794  result(i,j,0) = 0.0;
795  result(i,j,1) = 0.0;
796  result(i,j,2) = 0.0;
797  continue;
798  }
799 
800  auto m = cell_type.star(i,j);
801  auto U = velocity.star(i,j);
802 
803  // strain in units s-1
804  double u_x = 0, u_y = 0, v_x = 0, v_y = 0,
805  east = 1, west = 1, south = 1, north = 1;
806 
807  // Computes u_x using second-order centered finite differences written as
808  // weighted sums of first-order one-sided finite differences.
809  //
810  // Given the cell layout
811  // *----n----*
812  // | |
813  // | |
814  // w e
815  // | |
816  // | |
817  // *----s----*
818  // east == 0 if the east neighbor of the current cell is ice-free. In
819  // this case we use the left- (west-) sided difference.
820  //
821  // If both neighbors in the east-west (x) direction are ice-free the
822  // x-derivative is set to zero (see u_x, v_x initialization above).
823  //
824  // Similarly in y-direction.
825  if (ice_free(m.e)) {
826  east = 0;
827  }
828  if (ice_free(m.w)) {
829  west = 0;
830  }
831  if (ice_free(m.n)) {
832  north = 0;
833  }
834  if (ice_free(m.s)) {
835  south = 0;
836  }
837 
838  if (west + east > 0) {
839  u_x = 1.0 / (dx * (west + east)) * (west * (U.ij.u - U[West].u) + east * (U[East].u - U.ij.u));
840  v_x = 1.0 / (dx * (west + east)) * (west * (U.ij.v - U[West].v) + east * (U[East].v - U.ij.v));
841  }
842 
843  if (south + north > 0) {
844  u_y = 1.0 / (dy * (south + north)) * (south * (U.ij.u - U[South].u) + north * (U[North].u - U.ij.u));
845  v_y = 1.0 / (dy * (south + north)) * (south * (U.ij.v - U[South].v) + north * (U[North].v - U.ij.v));
846  }
847 
848  double nu = 0.0;
849  flow_law.effective_viscosity(hardness(i, j),
850  secondInvariant_2D(Vector2(u_x, v_x), Vector2(u_y, v_y)),
851  &nu, NULL);
852 
853  //get deviatoric stresses
854  result(i,j,0) = 2.0*nu*u_x;
855  result(i,j,1) = 2.0*nu*v_y;
856  result(i,j,2) = nu*(u_y+v_x);
857  }
858 }
859 
860 
861 } // end of namespace stressbalance
862 } // 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 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
std::shared_ptr< EnthalpyConverter > Ptr
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
IceModelVec2S latitude
Definition: Geometry.hh:43
IceModelVec2S longitude
Definition: Geometry.hh:44
IceModelVec2S cell_grounded_fraction
Definition: Geometry.hh:58
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S ice_area_specific_volume
Definition: Geometry.hh:54
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(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
stencils::Star< int > star(int i, int j) const
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
stencils::Star< T > star(int i, int j) const
Definition: IceModelVec2.hh:62
double * get_column(int i, int j)
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
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
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
unsigned int ndof() const
Returns the number of degrees of freedom per grid point.
Definition: iceModelVec.cc:132
void add(double alpha, const IceModelVec &x)
Result: v <- v + alpha * x. Calls VecAXPY.
Definition: iceModelVec.cc:234
void set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
Definition: iceModelVec.cc:175
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
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
void effective_viscosity(double hardness, double gamma, double *nu, double *dnu) const
Computes the regularized effective viscosity and its derivative with respect to the second invariant ...
Definition: FlowLaw.cc:163
void hardness_n(const double *enthalpy, const double *pressure, unsigned int n, double *result) const
Definition: FlowLaw.cc:122
double exponent() const
Definition: FlowLaw.cc:75
const IceModelVec2S * no_model_ice_thickness
const IceModelVec2S * water_column_pressure
const IceModelVec3 * enthalpy
void dump(const char *filename) const
const IceModelVec2S * fracture_density
const IceModelVec2S * basal_yield_stress
const IceModelVec2S * no_model_surface_elevation
const IceModelVec3 * age
const IceModelVec2V * bc_values
const IceModelVec2Int * bc_mask
const IceModelVec2S * basal_melt_rate
const IceModelVec2Int * no_model_mask
Shallow stress balance modifier (such as the non-sliding SIA).
Definition: SSB_Modifier.hh:39
Shallow stress balance (such as the SSA).
const IceModelVec3 & velocity_u() const
Get components of the the 3D velocity field.
StressBalance(IceGrid::ConstPtr g, std::shared_ptr< ShallowStressBalance > sb, std::shared_ptr< SSB_Modifier > ssb_mod)
std::shared_ptr< SSB_Modifier > m_modifier
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
const IceModelVec2S & basal_frictional_heating() const
Get the basal frictional heating.
const IceModelVec3 & volumetric_strain_heating() const
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
const IceModelVec2V & advective_velocity() const
Get the thickness-advective (SSA) 2D velocity.
const SSB_Modifier * modifier() const
Returns a pointer to a stress balance modifier implementation.
virtual void compute_volumetric_strain_heating(const Inputs &inputs)
Computes the volumetric strain heating using horizontal velocity.
std::shared_ptr< ShallowStressBalance > m_shallow_stress_balance
double max_diffusivity() const
Get the max diffusivity (for the adaptive time-stepping).
virtual void compute_vertical_velocity(const IceModelVec2CellType &mask, const IceModelVec3 &u, const IceModelVec3 &v, const IceModelVec2S *bmr, IceModelVec3 &result)
Compute vertical velocity using incompressibility of the ice.
void init()
Initialize the StressBalance object.
const ShallowStressBalance * shallow() const
Returns a pointer to a shallow stress balance solver implementation.
std::string stdout_report() const
Produce a report string for the standard output.
void update(const Inputs &inputs, bool full_update)
Update all the fields if (full_update), only update diffusive flux and max. diffusivity otherwise.
const IceModelVec2Stag & diffusive_flux() const
Get the diffusive (SIA) vertically-averaged flux on the staggered grid.
const IceModelVec3 & velocity_v() const
const IceModelVec3 & velocity_w() const
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
void append_time(const File &file, const Config &config, double time_seconds)
Prepare a file for output.
Definition: io_helpers.cc:256
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
Definition: io_helpers.cc:219
bool ice_free(int M)
Ice-free cell (grounded or ocean).
Definition: Mask.hh:57
void compute_2D_stresses(const rheology::FlowLaw &flow_law, const IceModelVec2V &velocity, const IceModelVec2S &hardness, const IceModelVec2CellType &cell_type, IceModelVec3 &result)
Compute 2D deviatoric stresses.
void compute_2D_principal_strain_rates(const IceModelVec2V &V, const IceModelVec2CellType &mask, IceModelVec3 &result)
Compute eigenvalues of the horizontal, vertically-integrated strain rate tensor.
static double D2(double u_x, double u_y, double u_z, double v_x, double v_y, double v_z)
static const double g
Definition: exactTestP.cc:39
IO_Backend string_to_backend(const std::string &backend)
Definition: File.cc:60
CFLData max_timestep_cfl_2d(const IceModelVec2S &ice_thickness, const IceModelVec2CellType &cell_type, const IceModelVec2V &velocity)
Compute the CFL constant associated to first-order upwinding for the sliding contribution to mass con...
CFLData max_timestep_cfl_3d(const IceModelVec2S &ice_thickness, const IceModelVec2CellType &cell_type, const IceModelVec3 &u3, const IceModelVec3 &v3, const IceModelVec3 &w3)
Compute the maximum velocities for time-stepping and reporting to user.
Definition: timestepping.cc:43
static const double k
Definition: exactTestP.cc:45
@ North
Definition: stencils.hh:24
@ East
Definition: stencils.hh:24
@ South
Definition: stencils.hh:24
@ West
Definition: stencils.hh:24
static double secondInvariant_2D(const Vector2 &U_x, const Vector2 &U_y)
Definition: FlowLaw.hh:44
@ PISM_READWRITE_MOVE
create a file for writing, move foo.nc to foo.nc~ if present
Definition: IO_Flags.hh:55
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49