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