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