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