PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
StressBalance_diagnostics.cc
Go to the documentation of this file.
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2022, 2023 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 <cmath> // std::pow, std::sqrt
20 #include <algorithm> // std::max
21 
22 
23 #include "pism/stressbalance/StressBalance_diagnostics.hh"
24 #include "pism/stressbalance/SSB_Modifier.hh"
25 #include "pism/stressbalance/ShallowStressBalance.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/pism_utilities.hh"
31 #include "pism/util/array/CellType.hh"
32 #include "pism/rheology/FlowLaw.hh"
33 #include "pism/rheology/FlowLawFactory.hh"
34 #include "pism/util/Context.hh"
35 
36 namespace pism {
37 namespace stressbalance {
38 
40  DiagnosticList result = {
41  {"bfrict", Diagnostic::Ptr(new PSB_bfrict(this))},
42  {"velbar_mag", Diagnostic::Ptr(new PSB_velbar_mag(this))},
43  {"flux", Diagnostic::Ptr(new PSB_flux(this))},
44  {"flux_mag", Diagnostic::Ptr(new PSB_flux_mag(this))},
45  {"velbase_mag", Diagnostic::Ptr(new PSB_velbase_mag(this))},
46  {"velsurf_mag", Diagnostic::Ptr(new PSB_velsurf_mag(this))},
47  {"uvel", Diagnostic::Ptr(new PSB_uvel(this))},
48  {"vvel", Diagnostic::Ptr(new PSB_vvel(this))},
49  {"strainheat", Diagnostic::Ptr(new PSB_strainheat(this))},
50  {"velbar", Diagnostic::Ptr(new PSB_velbar(this))},
51  {"velbase", Diagnostic::Ptr(new PSB_velbase(this))},
52  {"velsurf", Diagnostic::Ptr(new PSB_velsurf(this))},
53  {"wvel", Diagnostic::Ptr(new PSB_wvel(this))},
54  {"wvelbase", Diagnostic::Ptr(new PSB_wvelbase(this))},
55  {"wvelsurf", Diagnostic::Ptr(new PSB_wvelsurf(this))},
56  {"wvel_rel", Diagnostic::Ptr(new PSB_wvel_rel(this))},
57  {"strain_rates", Diagnostic::Ptr(new PSB_strain_rates(this))},
58  {"vonmises_stress", Diagnostic::Ptr(new PSB_vonmises_stress(this))},
59  {"deviatoric_stresses", Diagnostic::Ptr(new PSB_deviatoric_stresses(this))},
60  {"pressure", Diagnostic::Ptr(new PSB_pressure(this))},
61  {"tauxz", Diagnostic::Ptr(new PSB_tauxz(this))},
62  {"tauyz", Diagnostic::Ptr(new PSB_tauyz(this))}
63  };
64 
65  if (m_config->get_flag("output.ISMIP6")) {
66  result["velmean"] = Diagnostic::Ptr(new PSB_velbar(this));
67  result["zvelbase"] = Diagnostic::Ptr(new PSB_wvelbase(this));
68  result["zvelsurf"] = Diagnostic::Ptr(new PSB_wvelsurf(this));
69  }
70 
71  // add diagnostics from the shallow stress balance and the "modifier"
72  result = pism::combine(result, m_shallow_stress_balance->diagnostics());
73  result = pism::combine(result, m_modifier->diagnostics());
74 
75  return result;
76 }
77 
79  return pism::combine(m_shallow_stress_balance->ts_diagnostics(),
80  m_modifier->ts_diagnostics());
81 }
82 
84  : Diag<StressBalance>(m) {
85 
86  auto ismip6 = m_config->get_flag("output.ISMIP6");
87 
88  // set metadata:
89  m_vars = {{m_sys, ismip6 ? "xvelmean" : "ubar"},
90  {m_sys, ismip6 ? "yvelmean" : "vbar"}};
91 
92  m_vars[0]
93  .long_name("vertical mean of horizontal ice velocity in the X direction")
94  .standard_name("land_ice_vertical_mean_x_velocity")
95  .units("m s-1")
96  .output_units("m year-1");
97  m_vars[1]
98  .long_name("vertical mean of horizontal ice velocity in the Y direction")
99  .standard_name("land_ice_vertical_mean_y_velocity")
100  .units("m s-1")
101  .output_units("m year-1");
102 }
103 
104 std::shared_ptr<array::Array> PSB_velbar::compute_impl() const {
105  // get the thickness
106  const array::Scalar* thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
107 
108  // Compute the vertically-integrated horizontal ice flux:
109  auto result = array::cast<array::Vector>(PSB_flux(model).compute());
110 
111  // Override metadata set by the flux computation
112  result->metadata(0) = m_vars[0];
113  result->metadata(1) = m_vars[1];
114 
115  array::AccessScope list{thickness, result.get()};
116 
117  for (auto p = m_grid->points(); p; p.next()) {
118  const int i = p.i(), j = p.j();
119  double thk = (*thickness)(i,j);
120 
121  // Ice flux is masked already, but we need to check for division
122  // by zero anyway.
123  if (thk > 0.0) {
124  (*result)(i,j) /= thk;
125  } else {
126  (*result)(i,j) = 0.0;
127  }
128  }
129 
130  return result;
131 }
132 
134  : Diag<StressBalance>(m) {
135 
136  // set metadata:
137  m_vars = {{m_sys, "velbar_mag"}};
138 
139  m_vars[0]
140  .long_name("magnitude of vertically-integrated horizontal velocity of ice")
141  .units("m second-1")
142  .output_units("m year-1");
143 
144  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
145  m_vars[0]["valid_min"] = {0.0};
146 }
147 
148 std::shared_ptr<array::Array> PSB_velbar_mag::compute_impl() const {
149  auto result = allocate<array::Scalar>("velbar_mag");
150 
151  // compute vertically-averaged horizontal velocity:
152  auto velbar = array::cast<array::Vector>(PSB_velbar(model).compute());
153 
154  // compute its magnitude:
155  compute_magnitude(*velbar, *result);
156 
157  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
158 
159  // mask out ice-free areas:
160  apply_mask(*thickness, to_internal(m_fill_value), *result);
161 
162  return result;
163 }
164 
165 
167  : Diag<StressBalance>(m) {
168 
169  // set metadata:
170  m_vars = {{m_sys, "uflux"}, {m_sys, "vflux"}};
171  m_vars[0]
172  .long_name("Vertically integrated horizontal flux of ice in the X direction")
173  .units("m2 s-1")
174  .output_units("m2 year-1");
175  m_vars[1]
176  .long_name("Vertically integrated horizontal flux of ice in the Y direction")
177  .units("m2 s-1")
178  .output_units("m2 year-1");
179 }
180 
181 std::shared_ptr<array::Array> PSB_flux::compute_impl() const {
182  double H_threshold = m_config->get_number("geometry.ice_free_thickness_standard");
183 
184  auto result = allocate<array::Vector>("flux");
185 
186  // get the thickness
187  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
188 
189  const array::Array3D
190  &u3 = model->velocity_u(),
191  &v3 = model->velocity_v();
192 
193  array::AccessScope list{&u3, &v3, thickness, result.get()};
194 
195  const auto &z = m_grid->z();
196 
197  ParallelSection loop(m_grid->com);
198  try {
199  for (auto p = m_grid->points(); p; p.next()) {
200  const int i = p.i(), j = p.j();
201 
202  double H = (*thickness)(i,j);
203 
204  // an "ice-free" cell:
205  if (H < H_threshold) {
206  (*result)(i, j) = 0.0;
207  continue;
208  }
209 
210  // an icy cell:
211  {
212  auto u = u3.get_column(i, j);
213  auto v = v3.get_column(i, j);
214 
215  Vector2d Q(0.0, 0.0);
216 
217  // ks is "k just below the surface"
218  int ks = m_grid->kBelowHeight(H);
219 
220  if (ks > 0) {
221  Vector2d v0(u[0], v[0]);
222 
223  for (int k = 1; k <= ks; ++k) {
224  Vector2d v1(u[k], v[k]);
225 
226  // trapezoid rule
227  Q += (z[k] - z[k - 1]) * 0.5 * (v0 + v1);
228 
229  v0 = v1;
230  }
231  }
232 
233  // rectangle method to integrate over the last level
234  Q += (H - z[ks]) * Vector2d(u[ks], v[ks]);
235 
236  (*result)(i, j) = Q;
237  }
238  }
239  } catch (...) {
240  loop.failed();
241  }
242  loop.check();
243 
244  return result;
245 }
246 
247 
249  : Diag<StressBalance>(m) {
250 
251  // set metadata:
252  m_vars = { { m_sys, "flux_mag" } };
253  m_vars[0]
254  .long_name("magnitude of vertically-integrated horizontal flux of ice")
255  .units("m2 s-1")
256  .output_units("m2 year-1");
257  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
258  m_vars[0]["valid_min"] = {0.0};
259 }
260 
261 std::shared_ptr<array::Array> PSB_flux_mag::compute_impl() const {
262  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
263 
264  // Compute the vertically-averaged horizontal ice velocity:
265  auto result = array::cast<array::Scalar>(PSB_velbar_mag(model).compute());
266 
267  result->metadata() = m_vars[0];
268 
269  array::AccessScope list{thickness, result.get()};
270 
271  for (auto p = m_grid->points(); p; p.next()) {
272  const int i = p.i(), j = p.j();
273 
274  (*result)(i,j) *= (*thickness)(i,j);
275  }
276 
277  apply_mask(*thickness, to_internal(m_fill_value), *result);
278 
279  return result;
280 }
281 
283  : Diag<StressBalance>(m) {
284 
285  m_vars = { { m_sys, "velbase_mag" } };
286  m_vars[0]
287  .long_name("magnitude of horizontal velocity of ice at base of ice")
288  .units("m s-1")
289  .output_units("m year-1");
290  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
291  m_vars[0]["valid_min"] = {0.0};
292 }
293 
294 std::shared_ptr<array::Array> PSB_velbase_mag::compute_impl() const {
295  auto result = allocate<array::Scalar>("velbase_mag");
296 
297  compute_magnitude(*array::cast<array::Vector>(PSB_velbase(model).compute()), *result);
298 
299  double fill_value = to_internal(m_fill_value);
300 
301  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
302 
303  array::AccessScope list{&mask, result.get()};
304 
305  for (auto p = m_grid->points(); p; p.next()) {
306  const int i = p.i(), j = p.j();
307 
308  if (mask.ice_free(i, j)) {
309  (*result)(i, j) = fill_value;
310  }
311  }
312 
313  return result;
314 }
315 
317  : Diag<StressBalance>(m) {
318  m_vars = { { m_sys, "velsurf_mag" } };
319  m_vars[0]
320  .long_name("magnitude of horizontal velocity of ice at ice surface")
321  .units("m s-1")
322  .output_units("m year-1");
323  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
324  m_vars[0]["valid_min"] = {0.0};
325 }
326 
327 std::shared_ptr<array::Array> PSB_velsurf_mag::compute_impl() const {
328  double fill_value = to_internal(m_fill_value);
329 
330  auto result = allocate<array::Scalar>("velsurf_mag");
331 
332  compute_magnitude(*array::cast<array::Vector>(PSB_velsurf(model).compute()), *result);
333 
334  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
335 
336  array::AccessScope list{&mask, result.get()};
337 
338  for (auto p = m_grid->points(); p; p.next()) {
339  const int i = p.i(), j = p.j();
340 
341  if (mask.ice_free(i, j)) {
342  (*result)(i, j) = fill_value;
343  }
344  }
345 
346  return result;
347 }
348 
349 
351  : Diag<StressBalance>(m) {
352 
353  auto ismip6 = m_config->get_flag("output.ISMIP6");
354  m_vars = { { m_sys, ismip6 ? "xvelsurf" : "uvelsurf" },
355  { m_sys, ismip6 ? "yvelsurf" : "vvelsurf" } };
356  m_vars[0]
357  .long_name("x-component of the horizontal velocity of ice at ice surface")
358  .standard_name("land_ice_surface_x_velocity"); // InitMIP "standard" name
359  m_vars[1]
360  .long_name("y-component of the horizontal velocity of ice at ice surface")
361  .standard_name("land_ice_surface_y_velocity"); // InitMIP "standard" name
362 
363  auto large_number = to_internal(1e6);
364  for (auto &v : m_vars) {
365  v.units("m s-1").output_units("m year-1");
366  v["valid_range"] = {-large_number, large_number};
367  v["_FillValue"] = {to_internal(m_fill_value)};
368  }
369 }
370 
371 std::shared_ptr<array::Array> PSB_velsurf::compute_impl() const {
372  double fill_value = to_internal(m_fill_value);
373 
374  auto result = allocate<array::Vector>("surf");
375 
376  array::Scalar u_surf(m_grid, "u_surf");
377  array::Scalar v_surf(m_grid, "v_surf");
378 
379  const array::Array3D
380  &u3 = model->velocity_u(),
381  &v3 = model->velocity_v();
382 
383  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
384 
385  extract_surface(u3, *thickness, u_surf);
386  extract_surface(v3, *thickness, v_surf);
387 
388  const auto &cell_type = *m_grid->variables().get_2d_cell_type("mask");
389 
390  array::AccessScope list{ &cell_type, &u_surf, &v_surf, result.get() };
391 
392  for (auto p = m_grid->points(); p; p.next()) {
393  const int i = p.i(), j = p.j();
394 
395  if (cell_type.ice_free(i, j)) {
396  (*result)(i, j) = fill_value;
397  } else {
398  (*result)(i, j) = { u_surf(i, j), v_surf(i, j) };
399  }
400  }
401 
402  return result;
403 }
404 
406  m_vars = { { m_sys, "wvel", m_grid->z() } };
407  m_vars[0]
408  .long_name("vertical velocity of ice, relative to geoid")
409  .units("m s-1")
410  .output_units("m year-1");
411 
412  auto large_number = to_internal(1e6);
413  m_vars[0]["valid_range"] = { -large_number, large_number };
414 }
415 
416 std::shared_ptr<array::Array> PSB_wvel::compute(bool zero_above_ice) const {
417  std::shared_ptr<array::Array3D> result3(
418  new array::Array3D(m_grid, "wvel", array::WITHOUT_GHOSTS, m_grid->z()));
419  result3->metadata() = m_vars[0];
420 
421  const array::Scalar *bed, *uplift;
422  bed = m_grid->variables().get_2d_scalar("bedrock_altitude");
423  uplift = m_grid->variables().get_2d_scalar("tendency_of_bedrock_altitude");
424 
425  const array::Scalar &thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
426  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
427 
428  const array::Array3D &u3 = model->velocity_u(), &v3 = model->velocity_v(),
429  &w3 = model->velocity_w();
430 
431  array::AccessScope list{ &thickness, &mask, bed, &u3, &v3, &w3, uplift, result3.get() };
432 
433  const double ice_density = m_config->get_number("constants.ice.density"),
434  sea_water_density = m_config->get_number("constants.sea_water.density"),
435  R = ice_density / sea_water_density;
436 
437  ParallelSection loop(m_grid->com);
438  try {
439  for (auto p = m_grid->points(); p; p.next()) {
440  const int i = p.i(), j = p.j();
441 
442  const double *u = u3.get_column(i, j), *v = v3.get_column(i, j), *w = w3.get_column(i, j);
443  double *result = result3->get_column(i, j);
444 
445  int ks = m_grid->kBelowHeight(thickness(i, j));
446 
447  // in the ice:
448  if (mask.grounded(i, j)) {
449  const double bed_dx = diff_x_p(*bed, i, j), bed_dy = diff_y_p(*bed, i, j),
450  uplift_ij = (*uplift)(i, j);
451  for (int k = 0; k <= ks; k++) {
452  result[k] = w[k] + uplift_ij + u[k] * bed_dx + v[k] * bed_dy;
453  }
454 
455  } else { // floating
456  const double z_sl = R * thickness(i, j), w_sl = w3.interpolate(i, j, z_sl);
457 
458  for (int k = 0; k <= ks; k++) {
459  result[k] = w[k] - w_sl;
460  }
461  }
462 
463  // above the ice:
464  if (zero_above_ice) {
465  // set to zeros
466  for (unsigned int k = ks + 1; k < m_grid->Mz(); k++) {
467  result[k] = 0.0;
468  }
469  } else {
470  // extrapolate using the topmost value
471  for (unsigned int k = ks + 1; k < m_grid->Mz(); k++) {
472  result[k] = result[ks];
473  }
474  }
475  }
476  } catch (...) {
477  loop.failed();
478  }
479  loop.check();
480 
481  return result3;
482 }
483 
484 std::shared_ptr<array::Array> PSB_wvel::compute_impl() const {
485  return this->compute(true); // fill wvel above the ice with zeros
486 }
487 
489 
490  auto ismip6 = m_config->get_flag("output.ISMIP6");
491 
492  // set metadata:
493  m_vars = { { m_sys, ismip6 ? "zvelsurf" : "wvelsurf" } };
494 
495  m_vars[0]
496  .long_name("vertical velocity of ice at ice surface, relative to the geoid")
497  .standard_name("land_ice_surface_upward_velocity") // InitMIP "standard" name
498  .units("m s-1")
499  .output_units("m year-1");
500 
501  auto large_number = to_internal(1e6);
502  m_vars[0]["valid_range"] = { -large_number, large_number };
503  m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
504 }
505 
506 std::shared_ptr<array::Array> PSB_wvelsurf::compute_impl() const {
507  double fill_value = to_internal(m_fill_value);
508 
509  auto result = allocate<array::Scalar>("wvelsurf");
510 
511  // here "false" means "don't fill w3 above the ice surface with zeros"
512  auto w3 = array::cast<array::Array3D>(PSB_wvel(model).compute(false));
513 
514  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
515 
516  extract_surface(*w3, *thickness, *result);
517 
518  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
519 
520  array::AccessScope list{ &mask, result.get() };
521 
522  for (auto p = m_grid->points(); p; p.next()) {
523  const int i = p.i(), j = p.j();
524 
525  if (mask.ice_free(i, j)) {
526  (*result)(i, j) = fill_value;
527  }
528  }
529 
530  return result;
531 }
532 
534 
535  auto ismip6 = m_config->get_flag("output.ISMIP6");
536 
537  m_vars = { { m_sys, ismip6 ? "zvelbase" : "wvelbase" } };
538  m_vars[0]
539  .long_name("vertical velocity of ice at the base of ice, relative to the geoid")
540  .standard_name("land_ice_basal_upward_velocity") // InitMIP "standard" name
541  .units("m s-1")
542  .output_units("m year-1");
543 
544  auto large_number = to_internal(1e6);
545 
546  m_vars[0]["valid_range"] = { -large_number, large_number };
547  m_vars[0]["_FillValue"] = { to_internal(m_fill_value) };
548 }
549 
550 std::shared_ptr<array::Array> PSB_wvelbase::compute_impl() const {
551  double fill_value = to_internal(m_fill_value);
552 
553  auto result = allocate<array::Scalar>("wvelbase");
554 
555  // here "false" means "don't fill w3 above the ice surface with zeros"
556  auto w3 = array::cast<array::Array3D>(PSB_wvel(model).compute(false));
557 
558  extract_surface(*w3, 0.0, *result);
559 
560  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
561 
562  array::AccessScope list{ &mask, result.get() };
563 
564  for (auto p = m_grid->points(); p; p.next()) {
565  const int i = p.i(), j = p.j();
566 
567  if (mask.ice_free(i, j)) {
568  (*result)(i, j) = fill_value;
569  }
570  }
571 
572  return result;
573 }
574 
576 
577  auto ismip6 = m_config->get_flag("output.ISMIP6");
578 
579  // set metadata:
580  m_vars = { { m_sys, ismip6 ? "xvelbase" : "uvelbase" },
581  { m_sys, ismip6 ? "yvelbase" : "vvelbase" } };
582  m_vars[0]
583  .long_name("x-component of the horizontal velocity of ice at the base of ice")
584  .standard_name("land_ice_basal_x_velocity"); // InitMIP "standard" name
585  m_vars[1]
586  .long_name("y-component of the horizontal velocity of ice at the base of ice")
587  .standard_name("land_ice_basal_y_velocity"); // InitMIP "standard" name
588 
589  auto fill_value = to_internal(m_fill_value);
590  auto large_number = to_internal(1e6);
591 
592  for (auto &v : m_vars) {
593  v.units("m s-1").output_units("m year-1");
594  v["valid_range"] = { -large_number, large_number };
595  v["_FillValue"] = { fill_value };
596  }
597 }
598 
599 std::shared_ptr<array::Array> PSB_velbase::compute_impl() const {
600  double fill_value = to_internal(m_fill_value);
601 
602  auto result = allocate<array::Vector>("base");
603 
604  array::Scalar u_base(m_grid, "u_base");
605  array::Scalar v_base(m_grid, "v_base");
606 
607  const array::Array3D &u3 = model->velocity_u(), &v3 = model->velocity_v();
608 
609  extract_surface(u3, 0.0, u_base);
610  extract_surface(v3, 0.0, v_base);
611 
612  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
613 
614  array::AccessScope list{ &mask, &u_base, &v_base, result.get() };
615 
616  for (auto p = m_grid->points(); p; p.next()) {
617  const int i = p.i(), j = p.j();
618 
619  if (mask.ice_free(i, j)) {
620  (*result)(i, j) = fill_value;
621  } else {
622  (*result)(i, j) = { u_base(i, j), v_base(i, j) };
623  }
624  }
625 
626  return result;
627 }
628 
629 
631  m_vars = { { m_sys, "bfrict" } };
632  m_vars[0].long_name("basal frictional heating").units("W m-2");
633 }
634 
635 std::shared_ptr<array::Array> PSB_bfrict::compute_impl() const {
636 
637  auto result = allocate<array::Scalar>("bfrict");
638 
639  result->copy_from(model->basal_frictional_heating());
640 
641  return result;
642 }
643 
644 
646  m_vars = { { m_sys, "uvel", m_grid->z() } };
647  m_vars[0]
648  .long_name("horizontal velocity of ice in the X direction")
649  .standard_name("land_ice_x_velocity")
650  .units("m s-1")
651  .output_units("m year-1");
652 }
653 
654 /*!
655  * Copy F to result and set it to zero above the surface of the ice.
656  */
657 static void zero_above_ice(const array::Array3D &F, const array::Scalar &H,
658  array::Array3D &result) {
659 
660  array::AccessScope list{ &F, &H, &result };
661 
662  auto grid = result.grid();
663 
664  auto Mz = grid->Mz();
665 
666  ParallelSection loop(grid->com);
667  try {
668  for (auto p = grid->points(); p; p.next()) {
669  const int i = p.i(), j = p.j();
670 
671  int ks = grid->kBelowHeight(H(i, j));
672 
673  const double *F_ij = F.get_column(i, j);
674  double *F_out_ij = result.get_column(i, j);
675 
676  // in the ice:
677  for (int k = 0; k <= ks; k++) {
678  F_out_ij[k] = F_ij[k];
679  }
680  // above the ice:
681  for (unsigned int k = ks + 1; k < Mz; k++) {
682  F_out_ij[k] = 0.0;
683  }
684  }
685  } catch (...) {
686  loop.failed();
687  }
688  loop.check();
689 }
690 
691 std::shared_ptr<array::Array> PSB_uvel::compute_impl() const {
692 
693  std::shared_ptr<array::Array3D> result(
694  new array::Array3D(m_grid, "uvel", array::WITHOUT_GHOSTS, m_grid->z()));
695  result->metadata() = m_vars[0];
696 
697  zero_above_ice(model->velocity_u(), *m_grid->variables().get_2d_scalar("land_ice_thickness"),
698  *result);
699 
700  return result;
701 }
702 
704  m_vars = { { m_sys, "vvel", m_grid->z() } };
705  m_vars[0]
706  .long_name("horizontal velocity of ice in the Y direction")
707  .standard_name("land_ice_y_velocity")
708  .units("m s-1")
709  .output_units("m year-1");
710 }
711 
712 std::shared_ptr<array::Array> PSB_vvel::compute_impl() const {
713 
714  std::shared_ptr<array::Array3D> result(
715  new array::Array3D(m_grid, "vvel", array::WITHOUT_GHOSTS, m_grid->z()));
716  result->metadata() = m_vars[0];
717 
718  zero_above_ice(model->velocity_v(), *m_grid->variables().get_2d_scalar("land_ice_thickness"),
719  *result);
720 
721  return result;
722 }
723 
725  m_vars = { { m_sys, "wvel_rel", m_grid->z() } };
726  m_vars[0]
727  .long_name("vertical velocity of ice, relative to base of ice directly below")
728  .units("m s-1")
729  .output_units("m year-1");
730 }
731 
732 std::shared_ptr<array::Array> PSB_wvel_rel::compute_impl() const {
733 
734  std::shared_ptr<array::Array3D> result(
735  new array::Array3D(m_grid, "wvel_rel", array::WITHOUT_GHOSTS, m_grid->z()));
736  result->metadata() = m_vars[0];
737 
738  zero_above_ice(model->velocity_w(), *m_grid->variables().get_2d_scalar("land_ice_thickness"),
739  *result);
740 
741  return result;
742 }
743 
744 
746  m_vars = { { m_sys, "strainheat", m_grid->z() } };
747  m_vars[0]
748  .long_name("rate of strain heating in ice (dissipation heating)")
749  .units("W m-3")
750  .output_units("mW m-3");
751 }
752 
753 std::shared_ptr<array::Array> PSB_strainheat::compute_impl() const {
754  auto result =
755  std::make_shared<array::Array3D>(m_grid, "strainheat", array::WITHOUT_GHOSTS, m_grid->z());
756 
757  result->metadata() = m_vars[0];
758 
759  result->copy_from(model->volumetric_strain_heating());
760 
761  return result;
762 }
763 
765 
766  m_vars = { { m_sys, "eigen1" }, { m_sys, "eigen2" } };
767  m_vars[0]
768  .long_name("first eigenvalue of the horizontal, vertically-integrated strain rate tensor")
769  .units("s-1");
770  m_vars[1]
771  .long_name("second eigenvalue of the horizontal, vertically-integrated strain rate tensor")
772  .units("s-1");
773 }
774 
775 std::shared_ptr<array::Array> PSB_strain_rates::compute_impl() const {
776  auto velbar = array::cast<array::Vector>(PSB_velbar(model).compute());
777 
778  auto result = std::make_shared<array::Array2D<PrincipalStrainRates> >(m_grid, "strain_rates",
780  result->metadata(0) = m_vars[0];
781  result->metadata(1) = m_vars[1];
782 
783  array::Vector1 velbar_with_ghosts(m_grid, "velbar");
784 
785  // copy_from communicates ghosts
786  velbar_with_ghosts.copy_from(*velbar);
787 
788  array::CellType1 cell_type(m_grid, "cell_type");
789  {
790  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
791  cell_type.copy_from(mask);
792  }
793 
794  compute_2D_principal_strain_rates(velbar_with_ghosts, cell_type, *result);
795 
796  return result;
797 }
798 
800  m_vars = { { m_sys, "sigma_xx" }, { m_sys, "sigma_yy" }, { m_sys, "sigma_xy" } };
801 
802  m_vars[0].long_name("deviatoric stress in X direction").units("Pa");
803  m_vars[1].long_name("deviatoric stress in Y direction").units("Pa");
804  m_vars[2].long_name("deviatoric shear stress").units("Pa");
805 }
806 
807 std::shared_ptr<array::Array> PSB_deviatoric_stresses::compute_impl() const {
808 
809  auto result = std::make_shared<array::Array2D<stressbalance::DeviatoricStresses> >(
810  m_grid, "deviatoric_stresses", array::WITHOUT_GHOSTS);
811  result->metadata(0) = m_vars[0];
812  result->metadata(1) = m_vars[1];
813  result->metadata(2) = m_vars[2];
814 
815  const array::Array3D *enthalpy = m_grid->variables().get_3d_scalar("enthalpy");
816  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
817 
818  array::Scalar hardness(m_grid, "hardness");
819  array::Vector1 velocity(m_grid, "velocity");
820 
821  averaged_hardness_vec(*model->shallow()->flow_law(), *thickness, *enthalpy, hardness);
822 
823  // copy_from updates ghosts
824  velocity.copy_from(*array::cast<array::Vector>(PSB_velbar(model).compute()));
825 
826  array::CellType1 cell_type(m_grid, "cell_type");
827  {
828  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
829  cell_type.copy_from(mask);
830  }
831 
832  stressbalance::compute_2D_stresses(*model->shallow()->flow_law(), velocity, hardness, cell_type,
833  *result);
834 
835  return result;
836 }
837 
839  m_vars = { { m_sys, "pressure", m_grid->z() } };
840  m_vars[0].long_name("pressure in ice (hydrostatic)").units("Pa");
841 }
842 
843 std::shared_ptr<array::Array> PSB_pressure::compute_impl() const {
844 
845  std::shared_ptr<array::Array3D> result(
846  new array::Array3D(m_grid, "pressure", array::WITHOUT_GHOSTS, m_grid->z()));
847  result->metadata(0) = m_vars[0];
848 
849  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
850 
851  array::AccessScope list{ thickness, result.get() };
852 
853  const double rg = m_config->get_number("constants.ice.density") *
854  m_config->get_number("constants.standard_gravity");
855 
856  ParallelSection loop(m_grid->com);
857  try {
858  for (auto p = m_grid->points(); p; p.next()) {
859  const int i = p.i(), j = p.j();
860 
861  unsigned int ks = m_grid->kBelowHeight((*thickness)(i, j));
862  double *P_out_ij = result->get_column(i, j);
863  const double H = (*thickness)(i, j);
864  // within the ice:
865  for (unsigned int k = 0; k <= ks; ++k) {
866  P_out_ij[k] = rg * (H - m_grid->z(k)); // FIXME: add atmospheric pressure?
867  }
868  // above the ice:
869  for (unsigned int k = ks + 1; k < m_grid->Mz(); ++k) {
870  P_out_ij[k] = 0.0; // FIXME: use atmospheric pressure?
871  }
872  }
873  } catch (...) {
874  loop.failed();
875  }
876  loop.check();
877 
878  return result;
879 }
880 
881 
883  m_vars = {{m_sys, "tauxz", m_grid->z()}};
884  m_vars[0].long_name("shear stress xz component (in shallow ice approximation SIA)").units("Pa");
885 }
886 
887 
888 /*!
889  * The SIA-applicable shear stress component tauxz computed here is not used
890  * by the model. This implementation intentionally does not use the
891  * eta-transformation or special cases at ice margins.
892  * CODE DUPLICATION WITH PSB_tauyz
893  */
894 std::shared_ptr<array::Array> PSB_tauxz::compute_impl() const {
895 
896  std::shared_ptr<array::Array3D> result(
897  new array::Array3D(m_grid, "tauxz", array::WITHOUT_GHOSTS, m_grid->z()));
898  result->metadata() = m_vars[0];
899 
900  const array::Scalar *thickness, *surface;
901 
902  thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
903  surface = m_grid->variables().get_2d_scalar("surface_altitude");
904 
905  array::AccessScope list{ surface, thickness, result.get() };
906 
907  const double rg = m_config->get_number("constants.ice.density") *
908  m_config->get_number("constants.standard_gravity");
909 
910  ParallelSection loop(m_grid->com);
911  try {
912  for (auto p = m_grid->points(); p; p.next()) {
913  const int i = p.i(), j = p.j();
914 
915  unsigned int ks = m_grid->kBelowHeight((*thickness)(i, j));
916  double *tauxz_out_ij = result->get_column(i, j);
917  const double H = (*thickness)(i, j), dhdx = diff_x_p(*surface, i, j);
918 
919  // within the ice:
920  for (unsigned int k = 0; k <= ks; ++k) {
921  tauxz_out_ij[k] = -rg * (H - m_grid->z(k)) * dhdx;
922  }
923  // above the ice:
924  for (unsigned int k = ks + 1; k < m_grid->Mz(); ++k) {
925  tauxz_out_ij[k] = 0.0;
926  }
927  }
928  } catch (...) {
929  loop.failed();
930  }
931  loop.check();
932 
933  return result;
934 }
935 
936 
938  m_vars = { { m_sys, "tauyz", m_grid->z() } };
939  m_vars[0].long_name("shear stress yz component (in shallow ice approximation SIA)").units("Pa");
940 }
941 
942 
943 /*!
944  * The SIA-applicable shear stress component tauyz computed here is not used
945  * by the model. This implementation intentionally does not use the
946  * eta-transformation or special cases at ice margins.
947  * CODE DUPLICATION WITH PSB_tauxz
948  */
949 std::shared_ptr<array::Array> PSB_tauyz::compute_impl() const {
950 
951  std::shared_ptr<array::Array3D> result(
952  new array::Array3D(m_grid, "tauyz", array::WITHOUT_GHOSTS, m_grid->z()));
953  result->metadata(0) = m_vars[0];
954 
955  const array::Scalar *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
956  const array::Scalar *surface = m_grid->variables().get_2d_scalar("surface_altitude");
957 
958  array::AccessScope list{ surface, thickness, result.get() };
959 
960  const double rg = m_config->get_number("constants.ice.density") *
961  m_config->get_number("constants.standard_gravity");
962 
963  ParallelSection loop(m_grid->com);
964  try {
965  for (auto p = m_grid->points(); p; p.next()) {
966  const int i = p.i(), j = p.j();
967 
968  unsigned int ks = m_grid->kBelowHeight((*thickness)(i, j));
969  double *tauyz_out_ij = result->get_column(i, j);
970  const double H = (*thickness)(i, j), dhdy = diff_y_p(*surface, i, j);
971 
972  // within the ice:
973  for (unsigned int k = 0; k <= ks; ++k) {
974  tauyz_out_ij[k] = -rg * (H - m_grid->z(k)) * dhdy;
975  }
976  // above the ice:
977  for (unsigned int k = ks + 1; k < m_grid->Mz(); ++k) {
978  tauyz_out_ij[k] = 0.0;
979  }
980  }
981  } catch (...) {
982  loop.failed();
983  }
984  loop.check();
985 
986  return result;
987 }
988 
990  m_vars = { { m_sys, "vonmises_stress" } };
991  m_vars[0].long_name("tensile von Mises stress").units("Pascal");
992 }
993 
994 std::shared_ptr<array::Array> PSB_vonmises_stress::compute_impl() const {
995 
996  using std::max;
997  using std::sqrt;
998  using std::pow;
999 
1000  auto result = allocate<array::Scalar>("vonmises_stress");
1001 
1002  array::Scalar &vonmises_stress = *result;
1003 
1004  auto velbar = array::cast<array::Vector>(PSB_velbar(model).compute());
1005  array::Vector &velocity = *velbar;
1006 
1007  using StrainRates = array::Array2D<PrincipalStrainRates>;
1008  auto eigen12 = array::cast<StrainRates>(PSB_strain_rates(model).compute());
1009  const auto &strain_rates = *eigen12;
1010 
1011  const array::Scalar &ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
1012  const array::Array3D *enthalpy = m_grid->variables().get_3d_scalar("enthalpy");
1013  const auto &mask = *m_grid->variables().get_2d_cell_type("mask");
1014 
1015  std::shared_ptr<const rheology::FlowLaw> flow_law;
1016 
1017  if (m_config->get_flag("calving.vonmises_calving.use_custom_flow_law")) {
1018  EnthalpyConverter::Ptr EC = m_grid->ctx()->enthalpy_converter();
1019  rheology::FlowLawFactory factory("calving.vonmises_calving.", m_config, EC);
1020  flow_law = factory.create();
1021  } else {
1022  flow_law = model->shallow()->flow_law();
1023  }
1024 
1025  const double *z = &m_grid->z()[0];
1026 
1027  double glen_exponent = flow_law->exponent();
1028 
1029  array::AccessScope list{&vonmises_stress, &velocity, &strain_rates, &ice_thickness,
1030  enthalpy, &mask};
1031 
1032  for (auto pt = m_grid->points(); pt; pt.next()) {
1033  const int i = pt.i(), j = pt.j();
1034 
1035  if (mask.icy(i, j)) {
1036 
1037  const double H = ice_thickness(i, j);
1038  const unsigned int k = m_grid->kBelowHeight(H);
1039 
1040  const double
1041  *enthalpy_column = enthalpy->get_column(i, j),
1042  hardness = averaged_hardness(*flow_law, H, k, z, enthalpy_column),
1043  eigen1 = strain_rates(i, j).eigen1,
1044  eigen2 = strain_rates(i, j).eigen2;
1045 
1046  // [\ref Morlighem2016] equation 6
1047  const double effective_tensile_strain_rate = sqrt(0.5 * (pow(max(0.0, eigen1), 2) +
1048  pow(max(0.0, eigen2), 2)));
1049  // [\ref Morlighem2016] equation 7
1050  vonmises_stress(i, j) = sqrt(3.0) * hardness * pow(effective_tensile_strain_rate,
1051  1.0 / glen_exponent);
1052 
1053  } else { // end of "if (mask.icy(i, j))"
1054  vonmises_stress(i, j) = 0.0;
1055  }
1056  } // end of loop over grid points
1057 
1058  return result;
1059 }
1060 
1061 } // end of namespace stressbalance
1062 } // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const StressBalance * model
Definition: Diagnostic.hh:172
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:166
double m_fill_value
fill value (used often enough to justify storing it)
Definition: Diagnostic.hh:122
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:116
double to_internal(double x) const
Definition: Diagnostic.cc:59
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:120
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:65
std::shared_ptr< const Grid > m_grid
the grid
Definition: Diagnostic.hh:114
std::shared_ptr< array::Array > compute() const
Compute a diagnostic quantity and return a pointer to a newly-allocated Array.
Definition: Diagnostic.cc:131
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:118
std::shared_ptr< EnthalpyConverter > Ptr
void failed()
Indicates a failure of a parallel section.
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2d.hh:29
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
A storage vector combining related fields in a struct.
Definition: Array2D.hh:32
double * get_column(int i, int j)
Definition: Array3D.cc:120
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: Array3D.hh:33
std::shared_ptr< const Grid > grid() const
Definition: Array.cc:132
std::shared_ptr< FlowLaw > create()
virtual std::shared_ptr< array::Array > compute_impl() const
Computes basal frictional heating.
virtual std::shared_ptr< array::Array > compute_impl() const
Reports the vertically-integrated (2D) deviatoric stresses.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes flux_mag, the magnitude of vertically-integrated horizontal flux of ice.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes uflux and vflux, components of vertically-integrated horizontal flux of ice.
virtual std::shared_ptr< array::Array > compute_impl() const
Reports the pressure within the ice (3D).
virtual std::shared_ptr< array::Array > compute_impl() const
Reports the vertically-integrated (2D) principal strain rates.
virtual std::shared_ptr< array::Array > compute_impl() const
Reports the volumetric strain heating (3D).
virtual std::shared_ptr< array::Array > compute_impl() const
Reports the xz component of the shear stress within the ice (3D), according to the SIA formula.
virtual std::shared_ptr< array::Array > compute_impl() const
Reports the yz component of the shear stress within the ice (3D), according to the SIA formula.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes the x-component of the horizontal ice velocity.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes velbar_mag, the magnitude of vertically-integrated horizontal velocity of ice and masks out ...
virtual std::shared_ptr< array::Array > compute_impl() const
Computes the vertically-averaged ice velocity.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes velbase_mag, the magnitude of horizontal velocity of ice at base of ice and masks out ice-fr...
virtual std::shared_ptr< array::Array > compute_impl() const
Computes horizontal ice velocity at the base of ice.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes velsurf_mag, the magnitude of horizontal ice velocity at the surface.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes velsurf, the horizontal velocity of ice at ice surface.
std::shared_ptr< array::Array > compute_impl() const
virtual std::shared_ptr< array::Array > compute_impl() const
Computes the y-component of the horizontal ice velocity.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes vertical velocity of ice, relative to the bed directly below.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes vertical ice velocity (relative to the geoid).
virtual std::shared_ptr< array::Array > compute_impl() const
Computes wvelbase, the vertical velocity of ice at the base of ice.
virtual std::shared_ptr< array::Array > compute_impl() const
Computes wvelsurf, the vertical velocity of ice at ice surface.
std::shared_ptr< SSB_Modifier > m_modifier
virtual TSDiagnosticList ts_diagnostics_impl() const
std::shared_ptr< ShallowStressBalance > m_shallow_stress_balance
virtual DiagnosticList diagnostics_impl() const
The class defining PISM's interface to the shallow stress balance code.
void extract_surface(const Array3D &data, double z, Scalar &output)
Copies a horizontal slice at level z of an Array3D into output.
Definition: Array3D.cc:135
void apply_mask(const array::Scalar &M, double fill, array::Scalar &result)
Masks out all the areas where by setting them to fill.
Definition: Scalar.cc:79
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
double diff_x_p(const array::Scalar &array, int i, int j)
Returns the x-derivative at i,j approximated using centered finite differences. Respects grid periodi...
Definition: Scalar.cc:108
@ WITHOUT_GHOSTS
Definition: Array.hh:62
void compute_magnitude(const array::Vector &input, array::Scalar &result)
Definition: Scalar.cc:66
double diff_y_p(const array::Scalar &array, int i, int j)
Returns the y-derivative at i,j approximated using centered finite differences. Respects grid periodi...
Definition: Scalar.cc:129
void averaged_hardness_vec(const FlowLaw &ice, const array::Scalar &thickness, const array::Array3D &enthalpy, array::Scalar &result)
Definition: FlowLaw.cc:181
double averaged_hardness(const FlowLaw &ice, double thickness, unsigned int kbelowH, const double *zlevels, const double *enthalpy)
Computes vertical average of B(E, p) ice hardness, namely .
Definition: FlowLaw.cc:213
void compute_2D_principal_strain_rates(const array::Vector1 &V, const array::CellType1 &mask, array::Array2D< PrincipalStrainRates > &result)
Compute eigenvalues of the horizontal, vertically-integrated strain rate tensor.
static void zero_above_ice(const array::Array3D &F, const array::Scalar &H, array::Array3D &result)
void compute_2D_stresses(const rheology::FlowLaw &flow_law, const array::Vector1 &velocity, const array::Scalar &hardness, const array::CellType1 &cell_type, array::Array2D< DeviatoricStresses > &result)
Compute 2D deviatoric stresses.
static double F(double SL, double B, double H, double alpha)
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:343
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
static const double v0
Definition: exactTestP.cc:50
static const double k
Definition: exactTestP.cc:42
T combine(const T &a, const T &b)