PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
SSAFD.cc
Go to the documentation of this file.
1 // Copyright (C) 2004--2022 Constantine Khroulev, Ed Bueler and Jed Brown
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 <cassert>
20 #include <stdexcept>
21 
22 #include "SSAFD.hh"
23 #include "SSAFD_diagnostics.hh"
24 #include "pism/util/Mask.hh"
25 #include "pism/basalstrength/basal_resistance.hh"
26 #include "pism/util/pism_options.hh"
27 #include "pism/rheology/FlowLaw.hh"
28 #include "pism/util/Vars.hh"
29 #include "pism/util/IceGrid.hh"
30 #include "pism/util/Time.hh"
31 #include "pism/util/IceModelVec2CellType.hh"
32 #include "pism/stressbalance/StressBalance.hh"
33 #include "pism/geometry/Geometry.hh"
34 #include "pism/util/pism_utilities.hh"
35 #include "pism/util/petscwrappers/DM.hh"
36 #include "pism/util/petscwrappers/Vec.hh"
37 
38 namespace pism {
39 namespace stressbalance {
40 
41 SSAFD::KSPFailure::KSPFailure(const char* reason)
42  : RuntimeError(ErrorLocation(), std::string("SSAFD KSP (linear solver) failed: ") + reason){
43  // empty
44 }
45 
46 SSAFD::PicardFailure::PicardFailure(const std::string &message)
47  : RuntimeError(ErrorLocation(), "SSAFD Picard iterations failed: " + message) {
48  // empty
49 }
50 
52  return new SSAFD(g);
53 }
54 
55 /*!
56 Because the FD implementation of the SSA uses Picard iteration, a PETSc KSP
57 and Mat are used directly. In particular we set up \f$A\f$
58 (Mat m_A) and a \f$b\f$ (= Vec m_b) and iteratively solve
59 linear systems
60  \f[ A x = b \f]
61 where \f$x\f$ (= Vec SSAX). A PETSc SNES object is never created.
62  */
64  : SSA(grid),
65  m_hardness(grid, "hardness", WITHOUT_GHOSTS),
66  m_nuH(grid, "nuH", WITH_GHOSTS),
67  m_nuH_old(grid, "nuH_old", WITH_GHOSTS),
68  m_work(grid, "work_vector", WITH_GHOSTS,
69  2 /* stencil width */),
70  m_b(grid, "right_hand_side", WITHOUT_GHOSTS),
71  m_velocity_old(grid, "velocity_old", WITH_GHOSTS)
72 {
73 
74  m_velocity_old.set_attrs("internal",
75  "old SSA velocity field; used for re-trying with a different epsilon",
76  "m s-1", "m s-1", "", 0);
77 
78  auto units = pism::printf("Pa s^(1/%f)", m_flow_law->exponent());
79 
80  m_hardness.set_attrs("diagnostic",
81  "vertically-averaged ice hardness",
82  "1", "1",
83  "", 0);
85 
86  m_nuH.set_attrs("internal",
87  "ice thickness times effective viscosity",
88  "Pa s m", "Pa s m", "", 0);
89 
90  m_nuH_old.set_attrs("internal",
91  "ice thickness times effective viscosity (before an update)",
92  "Pa s m", "Pa s m", "", 0);
93 
94  m_work.set_attrs("internal",
95  "temporary storage used to compute nuH",
96  "", "", "", 0);
97 
98  m_scaling = 1.0e9; // comparable to typical beta for an ice stream;
99 
100  // The nuH viewer:
101  m_view_nuh = false;
102  m_nuh_viewer_size = 300;
103 
104  // PETSc objects and settings
105  {
106  PetscErrorCode ierr;
107  ierr = DMSetMatType(*m_da, MATAIJ);
108  PISM_CHK(ierr, "DMSetMatType");
109 
110  ierr = DMCreateMatrix(*m_da, m_A.rawptr());
111  PISM_CHK(ierr, "DMCreateMatrix");
112 
113  ierr = KSPCreate(m_grid->com, m_KSP.rawptr());
114  PISM_CHK(ierr, "KSPCreate");
115 
116  ierr = KSPSetOptionsPrefix(m_KSP, "ssafd_");
117  PISM_CHK(ierr, "KSPSetOptionsPrefix");
118 
119  // Use non-zero initial guess (i.e. SSA velocities from the last
120  // solve() call).
121  ierr = KSPSetInitialGuessNonzero(m_KSP, PETSC_TRUE);
122  PISM_CHK(ierr, "KSPSetInitialGuessNonzero");
123 
124  // Use the initial residual norm.
125  ierr = KSPConvergedDefaultSetUIRNorm(m_KSP);
126  PISM_CHK(ierr, "KSPConvergedDefaultSetUIRNorm");
127  }
128 }
129 
130 //! @note Uses `PetscErrorCode` *intentionally*.
132  PetscErrorCode ierr;
133  PC pc;
134 
135  ierr = KSPSetType(m_KSP, KSPGMRES);
136  PISM_CHK(ierr, "KSPSetType");
137 
138  ierr = KSPSetOperators(m_KSP, m_A, m_A);
139  PISM_CHK(ierr, "KSPSetOperators");
140 
141  // Get the PC from the KSP solver:
142  ierr = KSPGetPC(m_KSP, &pc);
143  PISM_CHK(ierr, "KSPGetPC");
144 
145  // Set the PC type:
146  ierr = PCSetType(pc, PCBJACOBI);
147  PISM_CHK(ierr, "PCSetType");
148 
149  // Process options:
150  ierr = KSPSetFromOptions(m_KSP);
151  PISM_CHK(ierr, "KSPSetFromOptions");
152 }
153 
154 //! @note Uses `PetscErrorCode` *intentionally*.
156  PetscErrorCode ierr;
157  PC pc, sub_pc;
158 
159  // Set parameters equivalent to
160  // -ksp_type gmres -ksp_norm_type unpreconditioned -ksp_pc_side right -pc_type asm -sub_pc_type lu
161 
162  ierr = KSPSetType(m_KSP, KSPGMRES);
163  PISM_CHK(ierr, "KSPSetType");
164 
165  ierr = KSPSetOperators(m_KSP, m_A, m_A);
166  PISM_CHK(ierr, "KSPSetOperators");
167 
168  // Switch to using the "unpreconditioned" norm.
169  ierr = KSPSetNormType(m_KSP, KSP_NORM_UNPRECONDITIONED);
170  PISM_CHK(ierr, "KSPSetNormType");
171 
172  // Switch to "right" preconditioning.
173  ierr = KSPSetPCSide(m_KSP, PC_RIGHT);
174  PISM_CHK(ierr, "KSPSetPCSide");
175 
176  // Get the PC from the KSP solver:
177  ierr = KSPGetPC(m_KSP, &pc);
178  PISM_CHK(ierr, "KSPGetPC");
179 
180  // Set the PC type:
181  ierr = PCSetType(pc, PCASM);
182  PISM_CHK(ierr, "PCSetType");
183 
184  // Set the sub-KSP object to "preonly"
185  KSP *sub_ksp;
186  ierr = PCSetUp(pc);
187  PISM_CHK(ierr, "PCSetUp");
188 
189  ierr = PCASMGetSubKSP(pc, NULL, NULL, &sub_ksp);
190  PISM_CHK(ierr, "PCASMGetSubKSP");
191 
192  ierr = KSPSetType(*sub_ksp, KSPPREONLY);
193  PISM_CHK(ierr, "KSPSetType");
194 
195  // Set the PC of the sub-KSP to "LU".
196  ierr = KSPGetPC(*sub_ksp, &sub_pc);
197  PISM_CHK(ierr, "KSPGetPC");
198 
199  ierr = PCSetType(sub_pc, PCLU);
200  PISM_CHK(ierr, "PCSetType");
201 
202  // Let the user override all this:
203  // Process options:
204  ierr = KSPSetFromOptions(m_KSP);
205  PISM_CHK(ierr, "KSPSetFromOptions");
206 }
207 
209  SSA::init_impl();
210 
211  m_log->message(2,
212  " [using the KSP-based finite difference implementation]\n");
213 
214  // options
215  options::Integer viewer_size("-ssa_nuh_viewer_size", "nuH viewer size",
217  m_nuh_viewer_size = viewer_size;
218  m_view_nuh = options::Bool("-ssa_view_nuh", "Enable the SSAFD nuH runtime viewer");
219 
220  if (m_config->get_flag("stress_balance.calving_front_stress_bc")) {
221  m_log->message(2,
222  " using PISM-PIK calving-front stress boundary condition ...\n");
223  }
224 
227 }
228 
229 //! \brief Computes the right-hand side ("rhs") of the linear problem for the
230 //! Picard iteration and finite-difference implementation of the SSA equations.
231 /*!
232 The right side of the SSA equations is just the driving stress term
233  \f[ - \rho g H \nabla h. \f]
234 The basal stress is put on the left side of the system. This method builds the
235 discrete approximation of the right side. For more about the discretization
236 of the SSA equations, see comments for assemble_matrix().
237 
238 The values of the driving stress on the i,j grid come from a call to
239 compute_driving_stress().
240 
241 In the case of Dirichlet boundary conditions, the entries on the right-hand side
242 come from known velocity values. The fields m_bc_values and m_bc_mask are used for
243 this.
244  */
245 void SSAFD::assemble_rhs(const Inputs &inputs) {
246  using mask::ice_free;
247  using mask::ice_free_land;
248  using mask::ice_free_ocean;
249 
250  const IceModelVec2S
251  &thickness = inputs.geometry->ice_thickness,
252  &bed = inputs.geometry->bed_elevation,
253  &surface = inputs.geometry->ice_surface_elevation,
254  &sea_level = inputs.geometry->sea_level_elevation,
255  *water_column_pressure = inputs.water_column_pressure;
256 
257  const double
258  dx = m_grid->dx(),
259  dy = m_grid->dy(),
260  standard_gravity = m_config->get_number("constants.standard_gravity"),
261  rho_ocean = m_config->get_number("constants.sea_water.density"),
262  rho_ice = m_config->get_number("constants.ice.density");
263 
264  // This constant is for debugging: simulations should not depend on the choice of
265  // velocity used in ice-free areas.
266  const Vector2 ice_free_velocity(0.0, 0.0);
267 
268  const bool
269  use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc"),
270  flow_line_mode = m_config->get_flag("stress_balance.ssa.fd.flow_line_mode");
271 
272  // FIXME: bedrock_boundary is a misleading name
273  bool bedrock_boundary = m_config->get_flag("stress_balance.ssa.dirichlet_bc");
274 
277  m_mask,
278  inputs.no_model_mask,
279  m_taud);
280 
282 
283  if (inputs.bc_values and inputs.bc_mask) {
284  list.add({inputs.bc_values, inputs.bc_mask});
285  }
286 
287  if (use_cfbc) {
288  list.add({&thickness, &bed, &surface, &m_mask, &sea_level});
289  }
290 
291  if (use_cfbc and water_column_pressure) {
292  list.add(*water_column_pressure);
293  }
294 
295  m_b.set(0.0);
296 
297  for (Points p(*m_grid); p; p.next()) {
298  const int i = p.i(), j = p.j();
299 
300  Vector2 taud = m_taud(i, j);
301 
302  if (flow_line_mode) {
303  // no cross-flow driving stress in the flow line mode
304  taud.v = 0.0;
305  }
306 
307  if (inputs.bc_values and inputs.bc_mask->as_int(i, j) == 1) {
308  m_b(i, j).u = m_scaling * (*inputs.bc_values)(i, j).u;
309  m_b(i, j).v = m_scaling * (*inputs.bc_values)(i, j).v;
310  continue;
311  }
312 
313  if (use_cfbc) {
314  double H_ij = thickness(i,j);
315 
316  auto M = m_mask.star(i, j);
317 
318  // Note: this sets velocities at both ice-free ocean and ice-free
319  // bedrock to zero. This means that we need to set boundary conditions
320  // at both ice/ice-free-ocean and ice/ice-free-bedrock interfaces below
321  // to be consistent.
322  if (ice_free(M.ij)) {
323  m_b(i, j) = m_scaling * ice_free_velocity;
324  continue;
325  }
326 
327  if (is_marginal(i, j, bedrock_boundary)) {
328  // weights at the west, east, south, and north cell faces
329  int W = 0, E = 0, S = 0, N = 0;
330  // direct neighbors
331  if (bedrock_boundary) {
332  if (ice_free_ocean(M.e))
333  E = 1;
334  if (ice_free_ocean(M.w))
335  W = 1;
336  if (ice_free_ocean(M.n))
337  N = 1;
338  if (ice_free_ocean(M.s))
339  S = 1;
340  } else {
341  if (ice_free(M.e))
342  E = 1;
343  if (ice_free(M.w))
344  W = 1;
345  if (ice_free(M.n))
346  N = 1;
347  if (ice_free(M.s))
348  S = 1;
349  }
350 
351  double
352  P_ice = 0.5 * rho_ice * standard_gravity * H_ij,
353  P_water = 0.0;
354 
355  if (water_column_pressure) {
356  P_water = (*water_column_pressure)(i, j);
357  } else {
358  P_water = pism::average_water_column_pressure(H_ij, bed(i, j), sea_level(i, j),
359  rho_ice, rho_ocean, standard_gravity);
360  }
361 
362  double delta_p = H_ij * (P_ice - P_water);
363 
364  if (grid_edge(*m_grid, i, j) and
365  not (flow_line_mode or mask::grounded(M.ij))) {
366  // In regional setups grounded ice may extend to the edge of the domain. This
367  // condition ensures that at a domain edge the ice behaves as if it extends past
368  // the edge without a change in geometry.
369  delta_p = 0.0;
370  }
371 
372  {
373  // fjord walls, nunataks, etc
374  //
375  // Override weights if we are at the margin and the grid cell on the other side
376  // of the interface is ice-free and above the level of the ice surface.
377  //
378  // This effectively sets the pressure difference at the corresponding interface
379  // to zero, which is exactly what we need.
380  auto b = bed.star(i, j);
381  double h = surface(i, j);
382 
383  if (ice_free(M.n) and b.n > h) {
384  N = 0;
385  }
386  if (ice_free(M.e) and b.e > h) {
387  E = 0;
388  }
389  if (ice_free(M.s) and b.s > h) {
390  S = 0;
391  }
392  if (ice_free(M.w) and b.w > h) {
393  W = 0;
394  }
395  }
396 
397  // Note that if the current cell is "marginal" but not a CFBC
398  // location, the following two lines are equaivalent to the "usual
399  // case" below.
400  //
401  // Note: signs below (+E, -W, etc) are explained by directions of outward
402  // normal vectors at corresponding cell faces.
403  m_b(i, j).u = taud.u + (E - W) * delta_p / dx;
404  m_b(i, j).v = taud.v + (N - S) * delta_p / dy;
405 
406  continue;
407  } // end of "if (is_marginal(i, j))"
408 
409  // If we reached this point, then CFBC are enabled, but we are in the
410  // interior of a sheet or shelf. See "usual case" below.
411 
412  } // end of "if (use_cfbc)"
413 
414  // usual case: use already computed driving stress
415  m_b(i, j) = taud;
416  }
417 }
418 
419 
420 //! \brief Assemble the left-hand side matrix for the KSP-based, Picard iteration,
421 //! and finite difference implementation of the SSA equations.
422 /*!
423 Recall the SSA equations are
424 \f{align*}
425  - 2 \left[\nu H \left(2 u_x + v_y\right)\right]_x
426  - \left[\nu H \left(u_y + v_x\right)\right]_y
427  - \tau_{(b)1} &= - \rho g H h_x, \\
428  - \left[\nu H \left(u_y + v_x\right)\right]_x
429  - 2 \left[\nu H \left(u_x + 2 v_y\right)\right]_y
430  - \tau_{(b)2} &= - \rho g H h_y,
431 \f}
432 where \f$u\f$ is the \f$x\f$-component of the velocity and \f$v\f$ is the
433 \f$y\f$-component of the velocity.
434 
435 The coefficient \f$\nu\f$ is the vertically-averaged effective viscosity.
436 (The product \f$\nu H\f$ is computed by compute_nuH_staggered().)
437 The Picard iteration idea is that, to solve the nonlinear equations in which
438 the effective viscosity depends on the velocity, we freeze the effective
439 viscosity using its value at the current estimate of the velocity and we solve
440 the linear equations which come from this viscosity. In abstract symbols, the
441 Picard iteration replaces the above nonlinear SSA equations by a sequence of
442 linear problems
443 
444 \f[ A(U^{(k)}) U^{(k+1)} = b \f]
445 
446 where \f$A(U)\f$ is the matrix from discretizing the SSA equations supposing
447 the viscosity is a function of known velocities \f$U\f$, and where \f$U^{(k)}\f$
448 denotes the \f$k\f$th iterate in the process. The current method assembles \f$A(U)\f$.
449 
450 For ice shelves \f$\tau_{(b)i} = 0\f$ [\ref MacAyealetal].
451 For ice streams with a basal till modelled as a plastic material,
452 \f$\tau_{(b)i} = - \tau_c u_i/|\mathbf{u}|\f$ where
453 \f$\mathbf{u} = (u,v)\f$, \f$|\mathbf{u}| = \left(u^2 + v^2\right)^{1/2}\f$,
454 where \f$\tau_c(t,x,y)\f$ is the yield stress of the till [\ref SchoofStream].
455 More generally, ice streams can be modeled with a pseudo-plastic basal till;
456 see IceModel::initBasalTillModel() and IceModel::updateYieldStressUsingBasalWater()
457 and reference [\ref BKAJS]. The pseudo-plastic till model includes all power law
458 sliding relations and the linearly-viscous model for sliding,
459 \f$\tau_{(b)i} = - \beta u_i\f$ where \f$\beta\f$ is the basal drag
460 (friction) parameter [\ref MacAyeal]. In any case, PISM assumes that the basal shear
461 stress can be factored this way, *even if the coefficient depends on the
462 velocity*, \f$\beta(u,v)\f$. Such factoring is possible even in the case of
463 (regularized) plastic till. This scalar coefficient \f$\beta\f$ is what is
464 returned by IceBasalResistancePlasticLaw::drag().
465 
466 Note that the basal shear stress appears on the \em left side of the
467 linear system we actually solve. We believe this is crucial, because
468 of its effect on the spectrum of the linear approximations of each
469 stage. The effect on spectrum is clearest in the linearly-viscous till
470 case but there seems to be an analogous effect in the plastic till case.
471 
472 This method assembles the matrix for the left side of the above SSA equations.
473 The numerical method is finite difference. Suppose we use difference notation
474 \f$\delta_{+x}f^{i,j} = f^{i+1,j}-f^{i,j}\f$,
475 \f$\delta_{-x}f^{i,j} = f^{i,j}-f^{i-1,j}\f$, and
476 \f$\Delta_{x}f^{i,j} = f^{i+1,j}-f^{i-1,j}\f$, and corresponding notation for
477 \f$y\f$ differences, and that we write \f$N = \nu H\f$ then the first of the
478 two "concrete" SSA equations above has this discretization:
479 \f{align*}
480 - &2 \frac{N^{i+\frac{1}{2},j}}{\Delta x} \left[2\frac{\delta_{+x}u^{i,j}}{\Delta x} + \frac{\Delta_{y} v^{i+1,j} + \Delta_{y} v^{i,j}}{4 \Delta y}\right] + 2 \frac{N^{i-\frac{1}{2},j}}{\Delta x} \left[2\frac{\delta_{-x}u^{i,j}}{\Delta x} + \frac{\Delta_y v^{i,j} + \Delta_y v^{i-1,j}}{4 \Delta y}\right] \\
481 &\qquad- \frac{N^{i,j+\frac{1}{2}}}{\Delta y} \left[\frac{\delta_{+y} u^{i,j}}{\Delta y} + \frac{\Delta_x v^{i,j+1} + \Delta_x v^{i,j}}{4 \Delta x}\right] + \frac{N^{i,j-\frac{1}{2}}}{\Delta y} \left[\frac{\delta_{-y}u^{i,j}}{\Delta y} + \frac{\Delta_x v^{i,j} + \Delta_x v^{i,j-1}}{4 \Delta x}\right] - \tau_{(b)1}^{i,j} = - \rho g H^{i,j} \frac{\Delta_x h^{i,j}}{2\Delta x}.
482 \f}
483 As a picture, see Figure \ref ssastencil.
484 
485 \image html ssastencil.png "\b ssastencil: Stencil for our finite difference discretization of the first of the two scalar SSA equations. Triangles show staggered grid points where N = nu * H is evaluated. Circles and squares show where u and v are approximated, respectively."
486 \anchor ssastencil
487 
488 It follows immediately that the matrix we assemble in the current method has
489 13 nonzeros entries per row because, for this first SSA equation, there are 5
490 grid values of \f$u\f$ and 8 grid values of \f$v\f$ used in this scheme. For
491 the second equation we also have 13 nonzeros per row.
492 
493 FIXME: document use of DAGetMatrix and MatStencil and MatSetValuesStencil
494 
495 */
496 void SSAFD::assemble_matrix(const Inputs &inputs,
497  bool include_basal_shear, Mat A) {
498  using mask::grounded_ice;
499  using mask::ice_free;
500  using mask::ice_free_land;
501  using mask::ice_free_ocean;
502  using mask::icy;
503 
504  PetscErrorCode ierr = 0;
505 
506  // shortcut:
507  const IceModelVec2V &vel = m_velocity;
508 
509  const IceModelVec2S
510  &thickness = inputs.geometry->ice_thickness,
511  &bed = inputs.geometry->bed_elevation,
512  &surface = inputs.geometry->ice_surface_elevation,
513  &grounded_fraction = inputs.geometry->cell_grounded_fraction,
514  &tauc = *inputs.basal_yield_stress;
515 
516  const double
517  dx = m_grid->dx(),
518  dy = m_grid->dy(),
519  beta_lateral_margin = m_config->get_number("basal_resistance.beta_lateral_margin"),
520  beta_ice_free_bedrock = m_config->get_number("basal_resistance.beta_ice_free_bedrock");
521 
522  const bool
523  // FIXME: bedrock_boundary is a misleading name
524  bedrock_boundary = m_config->get_flag("stress_balance.ssa.dirichlet_bc"),
525  flow_line_mode = m_config->get_flag("stress_balance.ssa.fd.flow_line_mode"),
526  use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc"),
527  replace_zero_diagonal_entries =
528  m_config->get_flag("stress_balance.ssa.fd.replace_zero_diagonal_entries");
529 
530  ierr = MatZeroEntries(A);
531  PISM_CHK(ierr, "MatZeroEntries");
532 
533  IceModelVec::AccessList list{&m_nuH, &tauc, &vel, &m_mask, &bed, &surface};
534 
535  if (inputs.bc_values && inputs.bc_mask) {
536  list.add(*inputs.bc_mask);
537  }
538 
539  const bool sub_gl = m_config->get_flag("geometry.grounded_cell_fraction");
540  if (sub_gl) {
541  list.add(grounded_fraction);
542  }
543 
544  // handles friction of the ice cell along ice-free bedrock margins when bedrock higher than ice
545  // surface (in simplified setups)
546  bool lateral_drag_enabled=m_config->get_flag("stress_balance.ssa.fd.lateral_drag.enabled");
547  if (lateral_drag_enabled) {
548  list.add({&thickness, &bed, &surface});
549  }
550  double lateral_drag_viscosity=m_config->get_number("stress_balance.ssa.fd.lateral_drag.viscosity");
551  double HminFrozen=0.0;
552 
553  /* matrix assembly loop */
554  ParallelSection loop(m_grid->com);
555  try {
556  for (Points p(*m_grid); p; p.next()) {
557  const int i = p.i(), j = p.j();
558 
559  // Handle the easy case: provided Dirichlet boundary conditions
560  if (inputs.bc_values && inputs.bc_mask && inputs.bc_mask->as_int(i,j) == 1) {
561  // set diagonal entry to one (scaled); RHS entry will be known velocity;
564  continue;
565  }
566 
567  /* Provide shorthand for the following staggered coefficients nu H:
568  * c_n
569  * c_w c_e
570  * c_s
571  */
572  // const
573  double c_w = m_nuH(i-1,j,0);
574  double c_e = m_nuH(i,j,0);
575  double c_s = m_nuH(i,j-1,1);
576  double c_n = m_nuH(i,j,1);
577 
578  if (lateral_drag_enabled) {
579  // if option is set, the viscosity at ice-bedrock boundary layer will
580  // be prescribed and is a temperature-independent free (user determined) parameter
581 
582  // direct neighbors
583  auto M = m_mask.star(i, j);
584  auto H = thickness.star(i, j);
585  auto b = bed.star(i, j);
586  double h = surface(i, j);
587 
588  if (H.ij > HminFrozen) {
589  if (b.w > h and ice_free_land(M.w)) {
590  c_w = lateral_drag_viscosity * 0.5 * (H.ij + H.w);
591  }
592  if (b.e > h and ice_free_land(M.e)) {
593  c_e = lateral_drag_viscosity * 0.5 * (H.ij + H.e);
594  }
595  if (b.n > h and ice_free_land(M.n)) {
596  c_n = lateral_drag_viscosity * 0.5 * (H.ij + H.n);
597  }
598  if (b.s > h and ice_free_land(M.s)) {
599  c_s = lateral_drag_viscosity * 0.5 * (H.ij + H.s);
600  }
601  }
602  }
603 
604  // We use DAGetMatrix to obtain the SSA matrix, which means that all 18
605  // non-zeros get allocated, even though we use only 13 (or 14). The
606  // remaining 5 (or 4) coefficients are zeros, but we set them anyway,
607  // because this makes the code easier to understand.
608  const int n_nonzeros = 18;
609  MatStencil row, col[n_nonzeros];
610 
611  // |-----+-----+---+-----+-----|
612  // | NW | NNW | N | NNE | NE |
613  // | WNW | | | | | ENE |
614  // | W |-----|-o-|-----| E |
615  // | WSW | | | | | ESE |
616  // | SW | SSW | S | SSE | SE |
617  // |-----+-----+---+-----+-----|
618  //
619  // We use compass rose notation for weights corresponding to interfaces between
620  // cells around the current one (i, j). Here N corresponds to the interface between
621  // the cell (i, j) and the one to the north of it.
622  int N = 1, E = 1, S = 1, W = 1;
623 
624  // Similarly, we use compass rose notation for weights used to switch between
625  // centered and one-sided finite differences. Here NNE is the interface between
626  // cells N and NE, ENE - between E and NE, etc.
627  int NNW = 1, NNE = 1, SSW = 1, SSE = 1;
628  int WNW = 1, ENE = 1, WSW = 1, ESE = 1;
629 
630  int M_ij = m_mask.as_int(i,j);
631 
632  if (use_cfbc) {
633  auto M = m_mask.box(i, j);
634 
635  // Note: this sets velocities at both ice-free ocean and ice-free
636  // bedrock to zero. This means that we need to set boundary conditions
637  // at both ice/ice-free-ocean and ice/ice-free-bedrock interfaces below
638  // to be consistent.
639  if (ice_free(M.ij)) {
642  continue;
643  }
644 
645  if (is_marginal(i, j, bedrock_boundary)) {
646  // If at least one of the following four conditions is "true", we're
647  // at a CFBC location.
648  if (bedrock_boundary) {
649 
650  if (ice_free_ocean(M.e))
651  E = 0;
652  if (ice_free_ocean(M.w))
653  W = 0;
654  if (ice_free_ocean(M.n))
655  N = 0;
656  if (ice_free_ocean(M.s))
657  S = 0;
658 
659  // decide whether to use centered or one-sided differences
660  if (ice_free_ocean(M.n) || ice_free_ocean(M.ne))
661  NNE = 0;
662  if (ice_free_ocean(M.e) || ice_free_ocean(M.ne))
663  ENE = 0;
664  if (ice_free_ocean(M.e) || ice_free_ocean(M.se))
665  ESE = 0;
666  if (ice_free_ocean(M.s) || ice_free_ocean(M.se))
667  SSE = 0;
668  if (ice_free_ocean(M.s) || ice_free_ocean(M.sw))
669  SSW = 0;
670  if (ice_free_ocean(M.w) || ice_free_ocean(M.sw))
671  WSW = 0;
672  if (ice_free_ocean(M.w) || ice_free_ocean(M.nw))
673  WNW = 0;
674  if (ice_free_ocean(M.n) || ice_free_ocean(M.nw))
675  NNW = 0;
676 
677  } else { // if (not bedrock_boundary)
678 
679  if (ice_free(M.e))
680  E = 0;
681  if (ice_free(M.w))
682  W = 0;
683  if (ice_free(M.n))
684  N = 0;
685  if (ice_free(M.s))
686  S = 0;
687 
688  // decide whether to use centered or one-sided differences
689  if (ice_free(M.n) || ice_free(M.ne))
690  NNE = 0;
691  if (ice_free(M.e) || ice_free(M.ne))
692  ENE = 0;
693  if (ice_free(M.e) || ice_free(M.se))
694  ESE = 0;
695  if (ice_free(M.s) || ice_free(M.se))
696  SSE = 0;
697  if (ice_free(M.s) || ice_free(M.sw))
698  SSW = 0;
699  if (ice_free(M.w) || ice_free(M.sw))
700  WSW = 0;
701  if (ice_free(M.w) || ice_free(M.nw))
702  WNW = 0;
703  if (ice_free(M.n) || ice_free(M.nw))
704  NNW = 0;
705 
706  } // end of the else clause following "if (bedrock_boundary)"
707  } // end of "if (is_marginal(i, j, bedrock_boundary))"
708  } // end of "if (use_cfbc)"
709 
710  /* begin Maxima-generated code */
711  const double dx2 = dx*dx, dy2 = dy*dy, d4 = 4*dx*dy, d2 = 2*dx*dy;
712 
713  /* Coefficients of the discretization of the first equation; u first, then v. */
714  double eq1[] = {
715  0, -c_n*N/dy2, 0,
716  -4*c_w*W/dx2, (c_n*N+c_s*S)/dy2+(4*c_e*E+4*c_w*W)/dx2, -4*c_e*E/dx2,
717  0, -c_s*S/dy2, 0,
718  c_w*W*WNW/d2+c_n*NNW*N/d4, (c_n*NNE*N-c_n*NNW*N)/d4+(c_w*W*N-c_e*E*N)/d2, -c_e*E*ENE/d2-c_n*NNE*N/d4,
719  (c_w*W*WSW-c_w*W*WNW)/d2+(c_n*W*N-c_s*W*S)/d4, (c_n*E*N-c_n*W*N-c_s*E*S+c_s*W*S)/d4+(c_e*E*N-c_w*W*N-c_e*E*S+c_w*W*S)/d2, (c_e*E*ENE-c_e*E*ESE)/d2+(c_s*E*S-c_n*E*N)/d4,
720  -c_w*W*WSW/d2-c_s*SSW*S/d4, (c_s*SSW*S-c_s*SSE*S)/d4+(c_e*E*S-c_w*W*S)/d2, c_e*E*ESE/d2+c_s*SSE*S/d4,
721  };
722 
723  /* Coefficients of the discretization of the second equation; u first, then v. */
724  double eq2[] = {
725  c_w*W*WNW/d4+c_n*NNW*N/d2, (c_n*NNE*N-c_n*NNW*N)/d2+(c_w*W*N-c_e*E*N)/d4, -c_e*E*ENE/d4-c_n*NNE*N/d2,
726  (c_w*W*WSW-c_w*W*WNW)/d4+(c_n*W*N-c_s*W*S)/d2, (c_n*E*N-c_n*W*N-c_s*E*S+c_s*W*S)/d2+(c_e*E*N-c_w*W*N-c_e*E*S+c_w*W*S)/d4, (c_e*E*ENE-c_e*E*ESE)/d4+(c_s*E*S-c_n*E*N)/d2,
727  -c_w*W*WSW/d4-c_s*SSW*S/d2, (c_s*SSW*S-c_s*SSE*S)/d2+(c_e*E*S-c_w*W*S)/d4, c_e*E*ESE/d4+c_s*SSE*S/d2,
728  0, -4*c_n*N/dy2, 0,
729  -c_w*W/dx2, (4*c_n*N+4*c_s*S)/dy2+(c_e*E+c_w*W)/dx2, -c_e*E/dx2,
730  0, -4*c_s*S/dy2, 0,
731  };
732 
733  /* i indices */
734  const int I[] = {
735  i-1, i, i+1,
736  i-1, i, i+1,
737  i-1, i, i+1,
738  i-1, i, i+1,
739  i-1, i, i+1,
740  i-1, i, i+1,
741  };
742 
743  /* j indices */
744  const int J[] = {
745  j+1, j+1, j+1,
746  j, j, j,
747  j-1, j-1, j-1,
748  j+1, j+1, j+1,
749  j, j, j,
750  j-1, j-1, j-1,
751  };
752 
753  /* component indices */
754  const int C[] = {
755  0, 0, 0,
756  0, 0, 0,
757  0, 0, 0,
758  1, 1, 1,
759  1, 1, 1,
760  1, 1, 1,
761  };
762  /* end Maxima-generated code */
763 
764  /* Dragging ice experiences friction at the bed determined by the
765  * IceBasalResistancePlasticLaw::drag() methods. These may be a plastic,
766  * pseudo-plastic, or linear friction law. Dragging is done implicitly
767  * (i.e. on left side of SSA eqns). */
768  double beta_u = 0.0, beta_v = 0.0;
769  if (include_basal_shear) {
770  double beta = 0.0;
771  if (grounded_ice(M_ij)) {
772  beta = m_basal_sliding_law->drag(tauc(i,j), vel(i,j).u, vel(i,j).v);
773  } else if (ice_free_land(M_ij)) {
774  // apply drag even in this case, to help with margins; note ice free
775  // areas already have a strength extension
776  beta = beta_ice_free_bedrock;
777  }
778  if (sub_gl) {
779  // reduce the basal drag at grid cells that are partially grounded:
780  if (icy(M_ij)) {
781  beta = grounded_fraction(i,j) * m_basal_sliding_law->drag(tauc(i,j), vel(i,j).u, vel(i,j).v);
782  }
783  }
784  beta_u = beta;
785  beta_v = beta;
786  }
787 
788  {
789  // Set very high basal drag *in the direction along the boundary* at locations
790  // bordering "fjord walls".
791 
792  auto M = m_mask.star(i, j);
793  auto b = bed.star(i, j);
794  double h = surface(i, j);
795 
796  if ((ice_free(M.n) and b.n > h) or (ice_free(M.s) and b.s > h)) {
797  beta_u += beta_lateral_margin;
798  }
799 
800  if ((ice_free(M.e) and b.e > h) or (ice_free(M.w) and b.w > h)) {
801  beta_v += beta_lateral_margin;
802  }
803  }
804 
805  // add beta to diagonal entries
806  eq1[4] += beta_u;
807  eq2[13] += beta_v;
808 
809  if (flow_line_mode) {
810  // set values corresponding to a trivial equation v = 0
811  for (int k = 0; k < n_nonzeros; ++k) {
812  eq2[k] = 0.0;
813  }
814  eq2[13] = m_scaling;
815  }
816 
817  // check diagonal entries:
818  const double eps = 1e-16;
819  if (fabs(eq1[4]) < eps) {
820  if (replace_zero_diagonal_entries) {
821  eq1[4] = beta_ice_free_bedrock;
822  } else {
823  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "first (X) equation in the SSAFD system:"
824  " zero diagonal entry at a regular (not Dirichlet B.C.)"
825  " location: i = %d, j = %d\n", i, j);
826  }
827  }
828  if (fabs(eq2[13]) < eps) {
829  if (replace_zero_diagonal_entries) {
830  eq2[13] = beta_ice_free_bedrock;
831  } else {
832  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "second (Y) equation in the SSAFD system:"
833  " zero diagonal entry at a regular (not Dirichlet B.C.)"
834  " location: i = %d, j = %d\n", i, j);
835  }
836  }
837 
838  row.i = i;
839  row.j = j;
840  for (int m = 0; m < n_nonzeros; m++) {
841  col[m].i = I[m];
842  col[m].j = J[m];
843  col[m].c = C[m];
844  }
845 
846  // set coefficients of the first equation:
847  row.c = 0;
848  ierr = MatSetValuesStencil(A, 1, &row, n_nonzeros, col, eq1, INSERT_VALUES);
849  PISM_CHK(ierr, "MatSetValuesStencil");
850 
851  // set coefficients of the second equation:
852  row.c = 1;
853  ierr = MatSetValuesStencil(A, 1, &row, n_nonzeros, col, eq2, INSERT_VALUES);
854  PISM_CHK(ierr, "MatSetValuesStencil");
855  } // i,j-loop
856  } catch (...) {
857  loop.failed();
858  }
859  loop.check();
860 
861  ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
862  PISM_CHK(ierr, "MatAssemblyBegin");
863 
864  ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
865  PISM_CHK(ierr, "MatAssemblyEnd");
866 #if (Pism_DEBUG==1)
867  ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
868  PISM_CHK(ierr, "MatSetOption");
869 #endif
870 }
871 
872 //! \brief Compute the vertically-averaged horizontal velocity from the shallow
873 //! shelf approximation.
874 /*!
875 This is the main procedure in the SSAFD. It manages the nonlinear solve process
876 and the Picard iteration.
877 
878 The outer loop (over index `k`) is the nonlinear iteration. In this loop the effective
879 viscosity is computed by compute_nuH_staggered() and then the linear system is
880 set up and solved.
881 
882 Specifically, we call the PETSc procedure KSPSolve() to solve the linear system.
883 Solving the linear system is also a loop, an iteration, but it occurs
884 inside KSPSolve(). The user has full control of the KSP solve through the PETSc
885 interface. The default choicess for KSP type `-ksp_type` and preconditioner type
886 `-pc_type` are GMRES(30) for the former and block Jacobi with ILU on the
887 blocks for the latter. The defaults usually work. These choices are important
888 but poorly understood. The eigenvalues of the linearized
889 SSA vary with ice sheet geometry and temperature in ways that are not well-studied.
890 Nonetheless these eigenvalues determine the convergence of
891 this (inner) linear iteration. A well-chosen preconditioner can put the eigenvalues
892 in the right place so that the KSP can converge quickly.
893 
894 The preconditioner will behave differently on different numbers of
895 processors. If the user wants the results of SSA calculations to be
896 independent of the number of processors, then `-pc_type none` could
897 be used, but performance will be poor.
898 
899 If you want to test different KSP methods, it may be helpful to see how many
900 iterations were necessary. Use `-ksp_monitor`.
901 Initial testing implies that CGS takes roughly half the iterations of
902 GMRES(30), but is not significantly faster because the iterations are each
903 roughly twice as slow. The outputs of PETSc options `-ksp_monitor_singular_value`,
904 `-ksp_compute_eigenvalues` and `-ksp_plot_eigenvalues -draw_pause N`
905 (the last holds plots for N seconds) may be useful to diagnose.
906 
907 The outer loop terminates when the effective viscosity times thickness
908 is no longer changing much, according to the tolerance set by the
909 option `-ssafd_picard_rtol`. The outer loop also terminates when a maximum
910 number of iterations is exceeded. We save the velocity from the last
911 time step in order to have a better estimate of the effective
912 viscosity than the u=v=0 result.
913 
914 In truth there is an "outer outer" loop (over index `l`). This attempts to
915 over-regularize the effective viscosity if the nonlinear iteration (the "outer"
916 loop over `k`) is not converging with the default regularization. The same
917 over-regularization is attempted if the KSP object reports that it has not
918 converged.
919 
920 (An alternative for recovery in the KSP diverged case, suggested by Jed, is to
921 revert to a direct linear solve, either for the whole domain (not scalable) or
922 on the subdomains. This recovery alternative requires a more nontrivial choice
923 but it may be worthwhile. Note the user can already do `-pc_type asm
924 -sub_pc_type lu` at the command line, forcing subdomain direct solves.)
925 
926 FIXME: update this doxygen comment
927 */
928 void SSAFD::solve(const Inputs &inputs) {
929 
930  // Store away old SSA velocity (it might be needed in case a solver
931  // fails).
933 
934  // These computations do not depend on the solution, so they need to
935  // be done once.
936  {
937  assemble_rhs(inputs);
938  compute_hardav_staggered(inputs);
939  }
940 
941  for (unsigned int k = 0; k < 3; ++k) {
942  try {
943  if (k == 0) {
944  // default strategy
945  picard_iteration(inputs, m_config->get_number("stress_balance.ssa.epsilon"), 1.0);
946 
947  break;
948  } else if (k == 1) {
949  // try underrelaxing the iteration
950  const double underrelax = m_config->get_number("stress_balance.ssa.fd.nuH_iter_failure_underrelaxation");
951  m_log->message(1,
952  " re-trying with effective viscosity under-relaxation (parameter = %.2f) ...\n",
953  underrelax);
954  picard_iteration(inputs, m_config->get_number("stress_balance.ssa.epsilon"), underrelax);
955 
956  break;
957  } else if (k == 2) {
958  // try over-regularization
960 
961  break;
962  } else {
963  // if we reached this, then all strategies above failed
964  write_system_petsc("all_strategies_failed");
965  throw RuntimeError(PISM_ERROR_LOCATION, "all SSAFD strategies failed");
966  }
967  } catch (PicardFailure &f) {
968  // proceed to the next strategy
969  }
970  }
971 
972  // Post-process velocities if the user asked for it:
973  if (m_config->get_flag("stress_balance.ssa.fd.brutal_sliding")) {
974  const double brutal_sliding_scaleFactor = m_config->get_number("stress_balance.ssa.fd.brutal_sliding_scale");
975  m_velocity.scale(brutal_sliding_scaleFactor);
976 
978  }
979 }
980 
981 void SSAFD::picard_iteration(const Inputs &inputs,
982  double nuH_regularization,
983  double nuH_iter_failure_underrelax) {
984 
986  // Give BJACOBI another shot if we haven't tried it enough yet
987 
988  try {
990  picard_manager(inputs, nuH_regularization,
991  nuH_iter_failure_underrelax);
992 
993  } catch (KSPFailure &f) {
994 
996 
997  m_log->message(1,
998  " re-trying using the Additive Schwarz preconditioner...\n");
999 
1000  pc_setup_asm();
1001 
1003 
1004  picard_manager(inputs, nuH_regularization,
1005  nuH_iter_failure_underrelax);
1006  }
1007 
1008  } else {
1009  // otherwise use ASM
1010  pc_setup_asm();
1011 
1012  picard_manager(inputs, nuH_regularization,
1013  nuH_iter_failure_underrelax);
1014  }
1015 }
1016 
1017 //! \brief Manages the Picard iteration loop.
1018 void SSAFD::picard_manager(const Inputs &inputs,
1019  double nuH_regularization,
1020  double nuH_iter_failure_underrelax) {
1021  PetscErrorCode ierr;
1022  double nuH_norm, nuH_norm_change;
1023  // ksp_iterations should be a PetscInt because it is used in the
1024  // KSPGetIterationNumber() call below
1025  PetscInt ksp_iterations, ksp_iterations_total = 0, outer_iterations;
1026  KSPConvergedReason reason;
1027 
1028  unsigned int max_iterations = static_cast<int>(m_config->get_number("stress_balance.ssa.fd.max_iterations"));
1029  double ssa_relative_tolerance = m_config->get_number("stress_balance.ssa.fd.relative_convergence");
1030  bool verbose = m_log->get_threshold() >= 2,
1031  very_verbose = m_log->get_threshold() > 2;
1032 
1033  // set the initial guess:
1035 
1036  m_stdout_ssa.clear();
1037 
1038  bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
1039 
1040  if (use_cfbc == true) {
1041  compute_nuH_staggered_cfbc(*inputs.geometry, nuH_regularization, m_nuH);
1042  } else {
1043  compute_nuH_staggered(*inputs.geometry, nuH_regularization, m_nuH);
1044  }
1046 
1047  // outer loop
1048  for (unsigned int k = 0; k < max_iterations; ++k) {
1049 
1050  if (very_verbose) {
1051  m_stdout_ssa += pism::printf(" %2d:", k);
1052  }
1053 
1054  // in preparation of measuring change of effective viscosity:
1056 
1057  // assemble (or re-assemble) matrix, which depends on updated viscosity
1058  assemble_matrix(inputs, true, m_A);
1059 
1060  if (very_verbose) {
1061 
1062  m_stdout_ssa += "A:";
1063  }
1064 
1065  // Call PETSc to solve linear system by iterative method; "inner iteration":
1066  ierr = KSPSetOperators(m_KSP, m_A, m_A);
1067  PISM_CHK(ierr, "KSPSetOperator");
1068 
1069  ierr = KSPSolve(m_KSP, m_b.vec(), m_velocity_global.vec());
1070  PISM_CHK(ierr, "KSPSolve");
1071 
1072  // Check if diverged; report to standard out about iteration
1073  ierr = KSPGetConvergedReason(m_KSP, &reason);
1074  PISM_CHK(ierr, "KSPGetConvergedReason");
1075 
1076  if (reason < 0) {
1077  // KSP diverged
1078  m_log->message(1,
1079  "PISM WARNING: KSPSolve() reports 'diverged'; reason = %d = '%s'\n",
1080  reason, KSPConvergedReasons[reason]);
1081 
1082  write_system_petsc("kspdivergederror");
1083 
1084  // Tell the caller that we failed. (The caller might try again,
1085  // though.)
1086  throw KSPFailure(KSPConvergedReasons[reason]);
1087  }
1088 
1089  // report on KSP success; the "inner" iteration is done
1090  ierr = KSPGetIterationNumber(m_KSP, &ksp_iterations);
1091  PISM_CHK(ierr, "KSPGetIterationNumber");
1092 
1093  ksp_iterations_total += ksp_iterations;
1094 
1095  if (very_verbose) {
1096  m_stdout_ssa += pism::printf("S:%d,%d: ", (int)ksp_iterations, reason);
1097  }
1098 
1099  // limit ice speed
1100  {
1101  auto max_speed = m_config->get_number("stress_balance.ssa.fd.max_speed", "m second-1");
1102  int high_speed_counter = 0;
1103 
1105 
1106  for (Points p(*m_grid); p; p.next()) {
1107  const int i = p.i(), j = p.j();
1108 
1109  auto speed = m_velocity_global(i, j).magnitude();
1110 
1111  if (speed > max_speed) {
1112  m_velocity_global(i, j) *= max_speed / speed;
1113  high_speed_counter += 1;
1114  }
1115  }
1116 
1117  high_speed_counter = GlobalSum(m_grid->com, high_speed_counter);
1118 
1119  if (high_speed_counter > 0) {
1120  m_log->message(2, " SSA speed was capped at %d locations\n", high_speed_counter);
1121  }
1122  }
1123 
1124  // Communicate so that we have stencil width for evaluation of effective
1125  // viscosity on next "outer" iteration (and geometry etc. if done):
1126  // Note that copy_from() updates ghosts of m_velocity.
1128 
1129  // update viscosity and check for viscosity convergence
1130  if (use_cfbc == true) {
1131  compute_nuH_staggered_cfbc(*inputs.geometry, nuH_regularization, m_nuH);
1132  } else {
1133  compute_nuH_staggered(*inputs.geometry, nuH_regularization, m_nuH);
1134  }
1135 
1136  if (nuH_iter_failure_underrelax != 1.0) {
1137  m_nuH.scale(nuH_iter_failure_underrelax);
1138  m_nuH.add(1.0 - nuH_iter_failure_underrelax, m_nuH_old);
1139  }
1140  compute_nuH_norm(nuH_norm, nuH_norm_change);
1141 
1143 
1144  if (very_verbose) {
1145  m_stdout_ssa += pism::printf("|nu|_2, |Delta nu|_2/|nu|_2 = %10.3e %10.3e\n",
1146  nuH_norm, nuH_norm_change/nuH_norm);
1147 
1148  // assume that high verbosity shows interest in immediate
1149  // feedback about SSA iterations
1150  m_log->message(2, m_stdout_ssa);
1151 
1152  m_stdout_ssa.clear();
1153  }
1154 
1155  outer_iterations = k + 1;
1156 
1157  if (nuH_norm == 0 || nuH_norm_change / nuH_norm < ssa_relative_tolerance) {
1158  goto done;
1159  }
1160 
1161  } // outer loop (k)
1162 
1163  // If we're here, it means that we exceeded max_iterations and still
1164  // failed.
1165 
1166  throw PicardFailure(pism::printf("effective viscosity not converged after %d iterations\n"
1167  "with nuH_regularization=%8.2e.",
1168  max_iterations, nuH_regularization));
1169 
1170  done:
1171 
1172  if (very_verbose) {
1173  auto tempstr = pism::printf("... =%5d outer iterations, ~%3.1f KSP iterations each\n",
1174  (int)outer_iterations,
1175  ((double) ksp_iterations_total) / outer_iterations);
1176  m_stdout_ssa += tempstr;
1177  } else if (verbose) {
1178  // at default verbosity, just record last nuH_norm_change and iterations
1179  auto tempstr = pism::printf("%5d outer iterations, ~%3.1f KSP iterations each\n",
1180  (int)outer_iterations,
1181  ((double) ksp_iterations_total) / outer_iterations);
1182 
1183  m_stdout_ssa += tempstr;
1184  }
1185 
1186  if (verbose) {
1187  m_stdout_ssa = " SSA: " + m_stdout_ssa;
1188  }
1189 }
1190 
1191 //! Old SSAFD recovery strategy: increase the SSA regularization parameter.
1193  // this has no units; epsilon goes up by this ratio when previous value failed
1194  const double DEFAULT_EPSILON_MULTIPLIER_SSA = 4.0;
1195  double nuH_regularization = m_config->get_number("stress_balance.ssa.epsilon");
1196  unsigned int k = 0, max_tries = 5;
1197 
1198  if (nuH_regularization <= 0.0) {
1199  throw PicardFailure("will not attempt over-regularization of nuH\n"
1200  "because nuH_regularization == 0.0.");
1201  }
1202 
1203  while (k < max_tries) {
1205  m_log->message(1,
1206  " re-trying with nuH_regularization multiplied by %8.2f...\n",
1207  DEFAULT_EPSILON_MULTIPLIER_SSA);
1208 
1209  nuH_regularization *= DEFAULT_EPSILON_MULTIPLIER_SSA;
1210 
1211  try {
1212  // 1.0 is the under-relaxation parameter
1213  picard_iteration(inputs, nuH_regularization, 1.0);
1214  // if this call succeeded, stop over-regularizing
1215  break;
1216  }
1217  catch (PicardFailure &f) {
1218  k += 1;
1219 
1220  if (k == max_tries) {
1221  throw PicardFailure("exceeded the max. number of nuH over-regularization attempts");
1222  }
1223  }
1224  }
1225 }
1226 
1227 //! \brief Compute the norm of nu H and the change in nu H.
1228 /*!
1229 Verification and PST experiments
1230 suggest that an \f$L^1\f$ criterion for convergence is best. For verification
1231 there seems to be little difference, presumably because the solutions are smooth
1232 and the norms are roughly equivalent on a subspace of smooth functions. For PST,
1233 the \f$L^1\f$ criterion gives faster runs with essentially the same results.
1234 Presumably that is because rapid (temporal and spatial) variation in
1235 \f$\nu H\f$ occurs at margins, occupying very few horizontal grid cells.
1236 For the significant (e.g.~in terms of flux) parts of the flow, it is o.k. to ignore
1237 a bit of bad behavior at these few places, and \f$L^1\f$ ignores it more than
1238 \f$L^2\f$ (much less \f$L^\infty\f$, which might not work at all).
1239  */
1240 void SSAFD::compute_nuH_norm(double &norm, double &norm_change) {
1241 
1242  const double area = m_grid->cell_area();
1243  const NormType MY_NORM = NORM_1;
1244 
1245  // Test for change in nu
1246  m_nuH_old.add(-1, m_nuH);
1247 
1248  std::vector<double>
1249  nuNorm = m_nuH.norm(MY_NORM),
1250  nuChange = m_nuH_old.norm(MY_NORM);
1251 
1252  nuChange[0] *= area;
1253  nuChange[1] *= area;
1254  nuNorm[0] *= area;
1255  nuNorm[1] *= area;
1256 
1257  norm_change = sqrt(PetscSqr(nuChange[0]) + PetscSqr(nuChange[1]));
1258  norm = sqrt(PetscSqr(nuNorm[0]) + PetscSqr(nuNorm[1]));
1259 }
1260 
1261 //! \brief Computes vertically-averaged ice hardness on the staggered grid.
1263  const IceModelVec2S
1264  &thickness = inputs.geometry->ice_thickness;
1265 
1266  const IceModelVec3 &enthalpy = *inputs.enthalpy;
1267 
1268  const double
1269  *E_ij = NULL,
1270  *E_offset = NULL;
1271 
1272  std::vector<double> E(m_grid->Mz());
1273 
1274  IceModelVec::AccessList list{&thickness, &enthalpy, &m_hardness, &m_mask};
1275 
1276  ParallelSection loop(m_grid->com);
1277  try {
1278  for (Points p(*m_grid); p; p.next()) {
1279  const int i = p.i(), j = p.j();
1280 
1281  E_ij = enthalpy.get_column(i,j);
1282  for (int o=0; o<2; o++) {
1283  const int oi = 1-o, oj=o;
1284  double H;
1285 
1286  if (m_mask.icy(i,j) && m_mask.icy(i+oi,j+oj)) {
1287  H = 0.5 * (thickness(i,j) + thickness(i+oi,j+oj));
1288  } else if (m_mask.icy(i,j)) {
1289  H = thickness(i,j);
1290  } else {
1291  H = thickness(i+oi,j+oj);
1292  }
1293 
1294  if (H == 0) {
1295  m_hardness(i,j,o) = -1e6; // an obviously impossible value
1296  continue;
1297  }
1298 
1299  E_offset = enthalpy.get_column(i+oi,j+oj);
1300  // build a column of enthalpy values a the current location:
1301  for (unsigned int k = 0; k < m_grid->Mz(); ++k) {
1302  E[k] = 0.5 * (E_ij[k] + E_offset[k]);
1303  }
1304 
1306  H, m_grid->kBelowHeight(H),
1307  &(m_grid->z()[0]), &E[0]);
1308  } // o
1309  } // loop over points
1310  } catch (...) {
1311  loop.failed();
1312  }
1313  loop.check();
1314 
1316 }
1317 
1318 
1319 /*! @brief Correct vertically-averaged hardness using a
1320  parameterization of the fracture-induced softening.
1321 
1322  See T. Albrecht, A. Levermann; Fracture-induced softening for
1323  large-scale ice dynamics; (2013), The Cryosphere Discussions 7;
1324  4501-4544; DOI:10.5194/tcd-7-4501-2013
1325 
1326  Note that this paper proposes an adjustment of the enhancement factor:
1327 
1328  \f[E_{\text{effective}} = E \cdot (1 - (1-\epsilon) \phi)^{-n}.\f]
1329 
1330  Let \f$E_{\text{effective}} = E\cdot C\f$, where \f$C\f$ is the
1331  factor defined by the formula above.
1332 
1333  Recall that the effective viscosity is defined by
1334 
1335  \f[\nu(D) = \frac12 B D^{(1-n)/(2n)}\f]
1336 
1337  and the viscosity form of the flow law is
1338 
1339  \f[\sigma'_{ij} = E_{\text{effective}}^{-\frac1n}2\nu(D) D_{ij}.\f]
1340 
1341  Then
1342 
1343  \f[\sigma'_{ij} = E_{\text{effective}}^{-\frac1n}BD^{(1-n)/(2n)}D_{ij}.\f]
1344 
1345  Using the fact that \f$E_{\text{effective}} = E\cdot C\f$, this can be rewritten as
1346 
1347  \f[\sigma'_{ij} = E^{-\frac1n} \left(C^{-\frac1n}B\right) D^{(1-n)/(2n)}D_{ij}.\f]
1348 
1349  So scaling the enhancement factor by \f$C\f$ is equivalent to scaling
1350  ice hardness \f$B\f$ by \f$C^{-\frac1n}\f$.
1351 */
1352 void SSAFD::fracture_induced_softening(const IceModelVec2S *fracture_density) {
1353  if (not fracture_density) {
1354  return;
1355  }
1356 
1357  const double
1358  epsilon = m_config->get_number("fracture_density.softening_lower_limit"),
1359  n_glen = m_flow_law->exponent();
1360 
1361  IceModelVec::AccessList list{&m_hardness, fracture_density};
1362 
1363  for (Points p(*m_grid); p; p.next()) {
1364  const int i = p.i(), j = p.j();
1365 
1366  for (int o=0; o<2; o++) {
1367  const int oi = 1-o, oj=o;
1368 
1369  const double
1370  // fracture density on the staggered grid:
1371  phi = 0.5 * ((*fracture_density)(i,j) + (*fracture_density)(i+oi,j+oj)),
1372  // the line below implements equation (6) in the paper
1373  softening = pow((1.0-(1.0-epsilon)*phi), -n_glen);
1374 
1375  m_hardness(i,j,o) *= pow(softening,-1.0/n_glen);
1376  }
1377  }
1378 }
1379 
1380 //! \brief Compute the product of ice thickness and effective viscosity (on the
1381 //! staggered grid).
1382 /*!
1383 In PISM the product \f$\nu H\f$ can be
1384  - constant, or
1385  - can be computed with a constant ice hardness \f$\bar B\f$ (temperature-independent)
1386  but with dependence of the viscosity on the strain rates, or
1387  - it can depend on the strain rates \e and have a vertically-averaged ice
1388  hardness depending on temperature or enthalpy.
1389 
1390 The flow law in ice stream and ice shelf regions must, for now, be a
1391 (temperature-dependent) Glen law. This is the only flow law we know how to
1392 invert to ``viscosity form''. More general forms like Goldsby-Kohlstedt are
1393 not yet inverted.
1394 
1395 The viscosity form of a Glen law is
1396  \f[ \nu(T^*,D) = \frac{1}{2} B(T^*) D^{(1/n)-1}\, D_{ij} \f]
1397 where
1398  \f[ D_{ij} = \frac{1}{2} \left(\frac{\partial U_i}{\partial x_j} +
1399  \frac{\partial U_j}{\partial x_i}\right) \f]
1400 is the strain rate tensor and \f$B\f$ is an ice hardness related to
1401 the ice softness \f$A(T^*)\f$ by
1402  \f[ B(T^*)=A(T^*)^{-1/n} \f]
1403 in the case of a temperature dependent Glen-type law. (Here \f$T^*\f$ is the
1404 pressure-adjusted temperature.)
1405 
1406 The effective viscosity is then
1407  \f[ \nu = \frac{\bar B}{2} \left[\left(\frac{\partial u}{\partial x}\right)^2 +
1408  \left(\frac{\partial v}{\partial y}\right)^2 +
1409  \frac{\partial u}{\partial x} \frac{\partial v}{\partial y} +
1410  \frac{1}{4} \left(\frac{\partial u}{\partial y}
1411  + \frac{\partial v}{\partial x}\right)^2
1412  \right]^{(1-n)/(2n)} \f]
1413 where in the temperature-dependent case
1414  \f[ \bar B = \frac{1}{H}\,\int_b^h B(T^*)\,dz\f]
1415 This integral is approximately computed by the trapezoid rule.
1416 
1417 The result of this procedure is \f$\nu H\f$, not just \f$\nu\f$, this it is
1418 a vertical integral, not average, of viscosity.
1419 
1420 The resulting effective viscosity times thickness is regularized by ensuring that
1421 its minimum is at least \f$\epsilon\f$. This regularization constant is an argument.
1422 
1423 In this implementation we set \f$\nu H\f$ to a constant anywhere the ice is
1424 thinner than a certain minimum. See SSAStrengthExtension and compare how this
1425 issue is handled when -cfbc is set.
1426 */
1428  double nuH_regularization,
1429  IceModelVec2Stag &result) {
1430 
1431  const IceModelVec2V &uv = m_velocity; // shortcut
1432 
1433  IceModelVec::AccessList list{&result, &uv, &m_hardness, &geometry.ice_thickness};
1434 
1435  double
1436  n_glen = m_flow_law->exponent(),
1437  nu_enhancement_scaling = 1.0 / pow(m_e_factor, 1.0 / n_glen);
1438 
1439  const double dx = m_grid->dx(), dy = m_grid->dy();
1440 
1441  for (int o=0; o<2; ++o) {
1442  const int oi = 1 - o, oj=o;
1443  for (Points p(*m_grid); p; p.next()) {
1444  const int i = p.i(), j = p.j();
1445 
1446  const double H = 0.5 * (geometry.ice_thickness(i,j) + geometry.ice_thickness(i+oi,j+oj));
1447 
1448  if (H < strength_extension->get_min_thickness()) {
1449  result(i,j,o) = strength_extension->get_notional_strength();
1450  continue;
1451  }
1452 
1453  double u_x, u_y, v_x, v_y;
1454  // Check the offset to determine how to differentiate velocity
1455  if (o == 0) {
1456  u_x = (uv(i+1,j).u - uv(i,j).u) / dx;
1457  u_y = (uv(i,j+1).u + uv(i+1,j+1).u - uv(i,j-1).u - uv(i+1,j-1).u) / (4*dy);
1458  v_x = (uv(i+1,j).v - uv(i,j).v) / dx;
1459  v_y = (uv(i,j+1).v + uv(i+1,j+1).v - uv(i,j-1).v - uv(i+1,j-1).v) / (4*dy);
1460  } else {
1461  u_x = (uv(i+1,j).u + uv(i+1,j+1).u - uv(i-1,j).u - uv(i-1,j+1).u) / (4*dx);
1462  u_y = (uv(i,j+1).u - uv(i,j).u) / dy;
1463  v_x = (uv(i+1,j).v + uv(i+1,j+1).v - uv(i-1,j).v - uv(i-1,j+1).v) / (4*dx);
1464  v_y = (uv(i,j+1).v - uv(i,j).v) / dy;
1465  }
1466 
1467  double nu = 0.0;
1468  m_flow_law->effective_viscosity(m_hardness(i,j,o),
1469  secondInvariant_2D({u_x, v_x}, {u_y, v_y}),
1470  &nu, NULL);
1471 
1472  result(i,j,o) = nu * H;
1473 
1474  // include the SSA enhancement factor; in most cases m_e_factor is 1
1475  result(i,j,o) *= nu_enhancement_scaling;
1476 
1477  // We ensure that nuH is bounded below by a positive constant.
1478  result(i,j,o) += nuH_regularization;
1479 
1480  } // i,j-loop
1481  } // o-loop
1482 
1483 
1484  // Some communication
1485  result.update_ghosts();
1486 }
1487 
1488 /**
1489  * @brief Compute the product of ice viscosity and thickness on the
1490  * staggered grid. Used when CFBC is enabled.
1491  *
1492  * @param[out] result nu*H product
1493  * @param[in] nuH_regularization regularization parameter (added to nu*H to keep it away from zero)
1494  *
1495  * @return 0 on success
1496  */
1498  double nuH_regularization,
1499  IceModelVec2Stag &result) {
1500 
1501  const IceModelVec2S &thickness = geometry.ice_thickness;
1502 
1503  const IceModelVec2V &uv = m_velocity; // shortcut
1504 
1505  double
1506  n_glen = m_flow_law->exponent(),
1507  nu_enhancement_scaling = 1.0 / pow(m_e_factor, 1.0 / n_glen);
1508 
1509  const double dx = m_grid->dx(), dy = m_grid->dy();
1510 
1512 
1513  assert(m_velocity.stencil_width() >= 2);
1514  assert(m_mask.stencil_width() >= 2);
1515  assert(m_work.stencil_width() >= 1);
1516 
1517  for (PointsWithGhosts p(*m_grid); p; p.next()) {
1518  const int i = p.i(), j = p.j();
1519 
1520  // x-derivative, i-offset
1521  {
1522  if (m_mask.icy(i,j) && m_mask.icy(i+1,j)) {
1523  m_work(i,j).u_x = (uv(i+1,j).u - uv(i,j).u) / dx; // u_x
1524  m_work(i,j).v_x = (uv(i+1,j).v - uv(i,j).v) / dx; // v_x
1525  m_work(i,j).w_i = 1.0;
1526  } else {
1527  m_work(i,j).u_x = 0.0;
1528  m_work(i,j).v_x = 0.0;
1529  m_work(i,j).w_i = 0.0;
1530  }
1531  }
1532 
1533  // y-derivative, j-offset
1534  {
1535  if (m_mask.icy(i,j) && m_mask.icy(i,j+1)) {
1536  m_work(i,j).u_y = (uv(i,j+1).u - uv(i,j).u) / dy; // u_y
1537  m_work(i,j).v_y = (uv(i,j+1).v - uv(i,j).v) / dy; // v_y
1538  m_work(i,j).w_j = 1.0;
1539  } else {
1540  m_work(i,j).u_y = 0.0;
1541  m_work(i,j).v_y = 0.0;
1542  m_work(i,j).w_j = 0.0;
1543  }
1544  }
1545  }
1546 
1547  list.add({&result, &m_hardness, &thickness});
1548 
1549  for (Points p(*m_grid); p; p.next()) {
1550  const int i = p.i(), j = p.j();
1551 
1552  double u_x, u_y, v_x, v_y, H, nu, W;
1553  // i-offset
1554  {
1555  if (m_mask.icy(i,j) && m_mask.icy(i+1,j)) {
1556  H = 0.5 * (thickness(i,j) + thickness(i+1,j));
1557  }
1558  else if (m_mask.icy(i,j)) {
1559  H = thickness(i,j);
1560  } else {
1561  H = thickness(i+1,j);
1562  }
1563 
1564  if (H >= strength_extension->get_min_thickness()) {
1565  u_x = m_work(i,j).u_x;
1566  v_x = m_work(i,j).v_x;
1567 
1568  W = m_work(i,j).w_j + m_work(i,j-1).w_j + m_work(i+1,j-1).w_j + m_work(i+1,j).w_j;
1569  if (W > 0) {
1570  u_y = 1.0/W * (m_work(i,j).u_y + m_work(i,j-1).u_y +
1571  m_work(i+1,j-1).u_y + m_work(i+1,j).u_y);
1572  v_y = 1.0/W * (m_work(i,j).v_y + m_work(i,j-1).v_y +
1573  m_work(i+1,j-1).v_y + m_work(i+1,j).v_y);
1574  } else {
1575  u_y = 0.0;
1576  v_y = 0.0;
1577  }
1578 
1579  m_flow_law->effective_viscosity(m_hardness(i,j,0),
1580  secondInvariant_2D(Vector2(u_x, v_x),
1581  Vector2(u_y, v_y)),
1582  &nu, NULL);
1583  result(i,j,0) = nu * H;
1584  } else {
1585  result(i,j,0) = strength_extension->get_notional_strength();
1586  }
1587  }
1588 
1589  // j-offset
1590  {
1591  if (m_mask.icy(i,j) && m_mask.icy(i,j+1)) {
1592  H = 0.5 * (thickness(i,j) + thickness(i,j+1));
1593  } else if (m_mask.icy(i,j)) {
1594  H = thickness(i,j);
1595  } else {
1596  H = thickness(i,j+1);
1597  }
1598 
1599  if (H >= strength_extension->get_min_thickness()) {
1600  u_y = m_work(i,j).u_y;
1601  v_y = m_work(i,j).v_y;
1602 
1603  W = m_work(i,j).w_i + m_work(i-1,j).w_i + m_work(i-1,j+1).w_i + m_work(i,j+1).w_i;
1604  if (W > 0.0) {
1605  u_x = 1.0/W * (m_work(i,j).u_x + m_work(i-1,j).u_x +
1606  m_work(i-1,j+1).u_x + m_work(i,j+1).u_x);
1607  v_x = 1.0/W * (m_work(i,j).v_x + m_work(i-1,j).v_x +
1608  m_work(i-1,j+1).v_x + m_work(i,j+1).v_x);
1609  } else {
1610  u_x = 0.0;
1611  v_x = 0.0;
1612  }
1613 
1614  m_flow_law->effective_viscosity(m_hardness(i,j,1),
1615  secondInvariant_2D(Vector2(u_x, v_x),
1616  Vector2(u_y, v_y)),
1617  &nu, NULL);
1618  result(i,j,1) = nu * H;
1619  } else {
1620  result(i,j,1) = strength_extension->get_notional_strength();
1621  }
1622  }
1623 
1624  // adjustments:
1625  for (unsigned int o = 0; o < 2; ++o) {
1626  // include the SSA enhancement factor; in most cases ssa_enhancement_factor is 1
1627  result(i,j,o) *= nu_enhancement_scaling;
1628 
1629  // We ensure that nuH is bounded below by a positive constant.
1630  result(i,j,o) += nuH_regularization;
1631  }
1632  }
1633 
1634  // Some communication
1635  result.update_ghosts();
1636 }
1637 
1638 //! Update the nuH viewer, which shows log10(nu H).
1640 
1641  if (not m_view_nuh) {
1642  return;
1643  }
1644 
1645  IceModelVec2S tmp(m_grid, "nuH", WITHOUT_GHOSTS);
1646  tmp.set_attrs("temporary",
1647  "log10 of (viscosity * thickness)",
1648  "Pa s m", "Pa s m", "", 0);
1649 
1650  IceModelVec::AccessList list{&m_nuH, &tmp};
1651 
1652  for (Points p(*m_grid); p; p.next()) {
1653  const int i = p.i(), j = p.j();
1654 
1655  double avg_nuH = 0.5 * (m_nuH(i,j,0) + m_nuH(i,j,1));
1656  if (avg_nuH > 1.0e14) {
1657  tmp(i,j) = log10(avg_nuH);
1658  } else {
1659  tmp(i,j) = 14.0;
1660  }
1661  }
1662 
1663  if (not m_nuh_viewer) {
1664  m_nuh_viewer.reset(new petsc::Viewer(m_grid->com, "nuH", m_nuh_viewer_size,
1665  m_grid->Lx(), m_grid->Ly()));
1666  }
1667 
1668  tmp.view({m_nuh_viewer});
1669 }
1670 
1671 void SSAFD::set_diagonal_matrix_entry(Mat A, int i, int j, int component,
1672  double value) {
1673  MatStencil row, col;
1674 
1675  row.i = i;
1676  row.j = j;
1677  row.c = component;
1678 
1679  col.i = i;
1680  col.j = j;
1681  col.c = component;
1682 
1683  PetscErrorCode ierr = MatSetValuesStencil(A, 1, &row, 1, &col, &value, INSERT_VALUES);
1684  PISM_CHK(ierr, "MatSetValuesStencil");
1685 }
1686 
1687 //! \brief Checks if a cell is near or at the ice front.
1688 /*!
1689  * You need to create IceModelVec::AccessList object and add mask to it.
1690  *
1691  * Note that a cell is a CFBC location of one of four direct neighbors is ice-free.
1692  *
1693  * If one of the diagonal neighbors is ice-free we don't use the CFBC, but we
1694  * do need to compute weights used in the SSA discretization (see
1695  * assemble_matrix()) to avoid differencing across interfaces between icy
1696  * and ice-free cells.
1697  *
1698  * This method ensures that checks in assemble_rhs() and assemble_matrix() are
1699  * consistent.
1700  */
1701 bool SSAFD::is_marginal(int i, int j, bool ssa_dirichlet_bc) {
1702 
1703  auto M = m_mask.box(i, j);
1704 
1705  using mask::ice_free;
1706  using mask::ice_free_ocean;
1707  using mask::icy;
1708 
1709  if (ssa_dirichlet_bc) {
1710  return icy(M.ij) &&
1711  (ice_free(M.e) || ice_free(M.w) || ice_free(M.n) || ice_free(M.s) ||
1712  ice_free(M.ne) || ice_free(M.se) || ice_free(M.nw) || ice_free(M.sw));
1713  } else {
1714  return icy(M.ij) &&
1715  (ice_free_ocean(M.e) || ice_free_ocean(M.w) ||
1716  ice_free_ocean(M.n) || ice_free_ocean(M.s) ||
1717  ice_free_ocean(M.ne) || ice_free_ocean(M.se) ||
1718  ice_free_ocean(M.nw) || ice_free_ocean(M.sw));
1719  }
1720 }
1721 
1722 void SSAFD::write_system_petsc(const std::string &namepart) {
1723  PetscErrorCode ierr;
1724 
1725  // write a file with a fixed filename; avoid zillions of files
1726  std::string filename = "SSAFD_" + namepart + ".petsc";
1727  m_log->message(1,
1728  " writing linear system to PETSc binary file %s ...\n", filename.c_str());
1729 
1730  petsc::Viewer viewer; // will be destroyed automatically
1731  ierr = PetscViewerBinaryOpen(m_grid->com, filename.c_str(), FILE_MODE_WRITE,
1732  viewer.rawptr());
1733  PISM_CHK(ierr, "PetscViewerBinaryOpen");
1734 
1735  ierr = MatView(m_A, viewer);
1736  PISM_CHK(ierr, "MatView");
1737 
1738  ierr = VecView(m_b.vec(), viewer);
1739  PISM_CHK(ierr, "VecView");
1740 }
1741 
1743  : Diag<SSAFD>(m) {
1744 
1745  // set metadata:
1746  m_vars = {SpatialVariableMetadata(m_sys, "nuH[0]"),
1747  SpatialVariableMetadata(m_sys, "nuH[1]")};
1748 
1749  set_attrs("ice thickness times effective viscosity, i-offset", "",
1750  "Pa s m", "kPa s m", 0);
1751  set_attrs("ice thickness times effective viscosity, j-offset", "",
1752  "Pa s m", "kPa s m", 1);
1753 }
1754 
1756 
1758  result->metadata(0) = m_vars[0];
1759  result->metadata(1) = m_vars[1];
1760 
1761  result->copy_from(model->integrated_viscosity());
1762 
1763  return result;
1764 }
1765 
1768 
1769  result["nuH"] = Diagnostic::Ptr(new SSAFD_nuH(this));
1770 
1771  return result;
1772 }
1773 
1775  return m_nuH;
1776 }
1777 
1778 
1779 } // end of namespace stressbalance
1780 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
const SSAFD * model
Definition: Diagnostic.hh:166
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:161
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:108
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:112
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:64
IceGrid::ConstPtr m_grid
the grid
Definition: Diagnostic.hh:106
void set_attrs(const std::string &long_name, const std::string &standard_name, const std::string &units, const std::string &glaciological_units, unsigned int N=0)
A method for setting common variable attributes.
Definition: Diagnostic.cc:132
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
IceModelVec2S cell_grounded_fraction
Definition: Geometry.hh:58
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S sea_level_elevation
Definition: Geometry.hh:50
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
virtual double drag(double tauc, double vx, double vy) const
Compute the drag coefficient for the basal shear stress.
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool icy(int i, int j) const
stencils::Box< int > box(int i, int j) const
int as_int(int i, int j) const
stencils::Star< int > star(int i, int j) const
void add(double alpha, const IceModelVec2S &x)
stencils::Star< double > star(int i, int j) const
std::shared_ptr< IceModelVec2Stag > Ptr
Definition: iceModelVec.hh:454
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
void copy_from(const IceModelVec2< T > &source)
Definition: IceModelVec2.hh:97
void add(double alpha, const IceModelVec2< T > &x)
Definition: IceModelVec2.hh:89
void copy_from(const IceModelVec3 &input)
double * get_column(int i, int j)
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:334
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void view(std::vector< std::shared_ptr< petsc::Viewer > > viewers) const
View a 2D vector field using existing PETSc viewers.
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
std::shared_ptr< IceModelVec > Ptr
Definition: iceModelVec.hh:206
petsc::Vec & vec() const
Definition: iceModelVec.cc:342
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
std::vector< double > norm(int n) const
Computes the norm of all the components of an IceModelVec.
Definition: iceModelVec.cc:769
void add(double alpha, const IceModelVec &x)
Result: v <- v + alpha * x. Calls VecAXPY.
Definition: iceModelVec.cc:234
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
void set_units_without_validation(const std::string &value)
double v
Definition: Vector2.hh:103
double u
Definition: Vector2.hh:103
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
const IceModelVec2S * water_column_pressure
const IceModelVec3 * enthalpy
const IceModelVec2S * fracture_density
const IceModelVec2S * basal_yield_stress
const IceModelVec2V * bc_values
const IceModelVec2Int * bc_mask
const IceModelVec2Int * no_model_mask
KSPFailure(const char *reason)
Definition: SSAFD.cc:41
PicardFailure(const std::string &message)
Definition: SSAFD.cc:46
virtual IceModelVec::Ptr compute_impl() const
Definition: SSAFD.cc:1755
SSAFD_nuH(const SSAFD *m)
Definition: SSAFD.cc:1742
Reports the nuH (viscosity times thickness) product on the staggered grid.
IceModelVec2V m_velocity_old
Definition: SSAFD.hh:114
virtual void compute_nuH_norm(double &norm, double &norm_change)
Compute the norm of nu H and the change in nu H.
Definition: SSAFD.cc:1240
IceModelVec2< Work > m_work
Definition: SSAFD.hh:107
virtual void init_impl()
Initialize a generic regular-grid SSA solver.
Definition: SSAFD.cc:208
virtual void picard_strategy_regularization(const Inputs &inputs)
Old SSAFD recovery strategy: increase the SSA regularization parameter.
Definition: SSAFD.cc:1192
std::shared_ptr< petsc::Viewer > m_nuh_viewer
Definition: SSAFD.hh:120
const IceModelVec2Stag & integrated_viscosity() const
Definition: SSAFD.cc:1774
virtual void update_nuH_viewers()
Update the nuH viewer, which shows log10(nu H).
Definition: SSAFD.cc:1639
virtual void compute_nuH_staggered_cfbc(const Geometry &geometry, double nuH_regularization, IceModelVec2Stag &result)
Compute the product of ice viscosity and thickness on the staggered grid. Used when CFBC is enabled.
Definition: SSAFD.cc:1497
IceModelVec2V m_b
Definition: SSAFD.hh:111
virtual void fracture_induced_softening(const IceModelVec2S *fracture_density)
Correct vertically-averaged hardness using a parameterization of the fracture-induced softening.
Definition: SSAFD.cc:1352
virtual void solve(const Inputs &inputs)
Compute the vertically-averaged horizontal velocity from the shallow shelf approximation.
Definition: SSAFD.cc:928
virtual void write_system_petsc(const std::string &namepart)
Definition: SSAFD.cc:1722
IceModelVec2Stag m_nuH_old
Definition: SSAFD.hh:91
IceModelVec2Stag m_hardness
Definition: SSAFD.hh:91
IceModelVec2Stag m_nuH
Definition: SSAFD.hh:91
virtual DiagnosticList diagnostics_impl() const
Definition: SSAFD.cc:1766
unsigned int m_default_pc_failure_max_count
Definition: SSAFD.hh:117
virtual void assemble_matrix(const Inputs &inputs, bool include_basal_shear, Mat A)
Assemble the left-hand side matrix for the KSP-based, Picard iteration, and finite difference impleme...
Definition: SSAFD.cc:496
virtual void pc_setup_asm()
Definition: SSAFD.cc:155
virtual bool is_marginal(int i, int j, bool ssa_dirichlet_bc)
Checks if a cell is near or at the ice front.
Definition: SSAFD.cc:1701
virtual void compute_nuH_staggered(const Geometry &geometry, double nuH_regularization, IceModelVec2Stag &result)
Compute the product of ice thickness and effective viscosity (on the staggered grid).
Definition: SSAFD.cc:1427
virtual void compute_hardav_staggered(const Inputs &inputs)
Computes vertically-averaged ice hardness on the staggered grid.
Definition: SSAFD.cc:1262
virtual void picard_iteration(const Inputs &inputs, double nuH_regularization, double nuH_iter_failure_underrelax)
Definition: SSAFD.cc:981
virtual void picard_manager(const Inputs &inputs, double nuH_regularization, double nuH_iter_failure_underrelax)
Manages the Picard iteration loop.
Definition: SSAFD.cc:1018
virtual void pc_setup_bjacobi()
Definition: SSAFD.cc:131
unsigned int m_default_pc_failure_count
Definition: SSAFD.hh:116
SSAFD(IceGrid::ConstPtr g)
Definition: SSAFD.cc:63
virtual void assemble_rhs(const Inputs &inputs)
Computes the right-hand side ("rhs") of the linear problem for the Picard iteration and finite-differ...
Definition: SSAFD.cc:245
void set_diagonal_matrix_entry(Mat A, int i, int j, int component, double value)
Definition: SSAFD.cc:1671
PISM's SSA solver: the finite difference implementation.
Definition: SSAFD.hh:34
double get_notional_strength() const
Returns strength = (viscosity times thickness).
Definition: SSA.cc:65
double get_min_thickness() const
Returns minimum thickness to trigger use of extension.
Definition: SSA.cc:70
virtual DiagnosticList diagnostics_impl() const
Definition: SSA.cc:367
IceModelVec2V m_taud
Definition: SSA.hh:141
virtual void compute_driving_stress(const IceModelVec2S &ice_thickness, const IceModelVec2S &surface_elevation, const IceModelVec2CellType &cell_type, const IceModelVec2Int *no_model_mask, IceModelVec2V &result) const
Compute the gravitational driving stress.
Definition: SSA.cc:237
SSAStrengthExtension * strength_extension
Definition: SSA.hh:115
std::string m_stdout_ssa
Definition: SSA.hh:143
std::shared_ptr< petsc::DM > m_da
Definition: SSA.hh:146
IceModelVec2V m_velocity_global
Definition: SSA.hh:147
IceModelVec2CellType m_mask
Definition: SSA.hh:140
virtual void init_impl()
Initialize a generic regular-grid SSA solver.
Definition: SSA.cc:126
PISM's SSA solver.
Definition: SSA.hh:110
std::shared_ptr< rheology::FlowLaw > m_flow_law
IceBasalResistancePlasticLaw * m_basal_sliding_law
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
const double rho_ice
Definition: exactTestK.c:31
const double phi
Definition: exactTestK.c:41
static Vector3 row(const double A[3][3], size_t k)
Definition: Element.cc:46
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:47
bool grounded_ice(int M)
Definition: Mask.hh:50
bool ice_free_land(int M)
Definition: Mask.hh:63
bool ice_free_ocean(int M)
Definition: Mask.hh:60
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:43
bool ice_free(int M)
Ice-free cell (grounded or ocean).
Definition: Mask.hh:57
bool Bool(const std::string &option, const std::string &description)
Definition: options.cc:240
double averaged_hardness(const FlowLaw &ice, double thickness, int kbelowH, const double *zlevels, const double *enthalpy)
Computes vertical average of B(E, p) ice hardness, namely .
Definition: FlowLaw.cc:214
SSA * SSAFDFactory(IceGrid::ConstPtr g)
Constructs a new SSAFD.
Definition: SSAFD.cc:51
double average_water_column_pressure(double ice_thickness, double bed, double floatation_level, double rho_ice, double rho_water, double g)
static const double g
Definition: exactTestP.cc:39
std::string printf(const char *format,...)
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
static const double k
Definition: exactTestP.cc:45
static double secondInvariant_2D(const Vector2 &U_x, const Vector2 &U_y)
Definition: FlowLaw.hh:44
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
bool grid_edge(const IceGrid &grid, int i, int j)
Definition: IceGrid.hh:350
const int J[]
Definition: ssafd_code.cc:34
double eq1[]
Definition: ssafd_code.cc:4
const int I[]
Definition: ssafd_code.cc:24
const double d2
Definition: ssafd_code.cc:1
const int C[]
Definition: ssafd_code.cc:44
const double d4
Definition: ssafd_code.cc:1
const double dy2
Definition: ssafd_code.cc:1
double eq2[]
Definition: ssafd_code.cc:14
const double dx2
Definition: ssafd_code.cc:1
static double S(unsigned n)
Definition: test_cube.c:58