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