PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Blatter.cc
Go to the documentation of this file.
1 /* Copyright (C) 2020, 2021 PISM Authors
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 
20 #include <cassert> // assert
21 #include <cmath> // std::pow, std::fabs, std::log10
22 #include <algorithm> // std::max
23 #include <cstring> // memset
24 
25 #include "Blatter.hh"
26 #include "pism/util/fem/FEM.hh"
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/Vector2.hh"
29 
30 #include "util/DataAccess.hh"
31 #include "util/grid_hierarchy.hh"
32 #include "pism/util/node_types.hh"
33 
34 #include "pism/rheology/FlowLawFactory.hh"
35 
36 #include "pism/stressbalance/StressBalance.hh"
37 #include "pism/geometry/Geometry.hh"
38 #include "pism/util/pism_options.hh"
39 
40 namespace pism {
41 namespace stressbalance {
42 
43 /*!
44  * Compute node type using domain thickness and the thickness threshold `min_thickness`.
45  *
46  * An element contains ice if ice thickness at all its nodes equal or exceeds the
47  * `min_thickness` threshold.
48  *
49  * A node is *interior* if all four elements it belongs to contain ice.
50  *
51  * A node is *exterior* if it belongs to zero icy elements.
52  *
53  * A node that is neither interior nor exterior is a *boundary* node.
54  */
55 void Blatter::compute_node_type(double min_thickness) {
56 
57  IceModelVec2S node_type(m_grid, "node_type", WITH_GHOSTS);
58  node_type.set(0.0);
59 
60  DMDALocalInfo info;
61  int ierr = DMDAGetLocalInfo(m_da, &info); PISM_CHK(ierr, "DMDAGetLocalInfo");
62  info = grid_transpose(info);
63 
64  // Note that dx, dy, and quadrature don't matter here.
65  fem::Q1Element2 E(info, 1.0, 1.0, fem::Q1Quadrature1());
66 
68 
69  IceModelVec::AccessList l{&node_type, &m_parameters};
70 
71  // Loop over all the elements with at least one owned node and compute the number of icy
72  // elements each node belongs to.
73  for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
74  for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
75  E.reset(i, j);
76 
77  E.nodal_values(m_parameters.array(), p);
78 
79  // An element is "interior" (contains ice) if all of its nodes have thickness above
80  // the threshold
81  bool interior = true;
82  for (int k = 0; k < fem::q1::n_chi; ++k) {
83  if (p[k].thickness < min_thickness) {
84  interior = false;
85  break;
86  }
87  }
88 
89  for (int k = 0; k < fem::q1::n_chi; ++k) {
90  int ii, jj;
91  E.local_to_global(k, ii, jj);
92  node_type(ii, jj) += static_cast<double>(interior);
93  }
94  }
95  }
96 
97  node_type.update_ghosts();
98 
99  // Loop over all the owned nodes and turn the number of "icy" elements this node belongs
100  // to into node type.
101  for (int j = info.ys; j < info.ys + info.ym; j++) {
102  for (int i = info.xs; i < info.xs + info.xm; i++) {
103 
104  switch ((int)node_type(i, j)) {
105  case 4:
106  m_parameters(i, j).node_type = NODE_INTERIOR;
107  break;
108  case 0:
109  m_parameters(i, j).node_type = NODE_EXTERIOR;
110  break;
111  default:
112  m_parameters(i, j).node_type = NODE_BOUNDARY;
113  }
114  }
115  }
116 }
117 
118 /*!
119  * Returns true if a node is in the Dirichlet part of the boundary, false otherwise.
120  *
121  * Used by verification tests.
122  */
123 bool Blatter::dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex& I) {
124  (void) info;
125  (void) I;
126  return false;
127 }
128 
129 /*! Dirichlet BC
130 */
131 Vector2 Blatter::u_bc(double x, double y, double z) const {
132  (void) x;
133  (void) y;
134  (void) z;
135  return {0.0, 0.0};
136 }
137 
138 /*!
139  * Return true if an element does not contain ice, i.e. is a part of the "exterior" of the
140  * ice mass.
141  *
142  * @param[in] node_type node type at the nodes of an element (an array of 8 integers; only
143  * 4 are used)
144  */
145 bool Blatter::exterior_element(const int *node_type) {
146  // number of nodes per map-plane cell
147  int N = 4;
148  for (int n = 0; n < N; ++n) {
149  if (node_type[n] == NODE_EXTERIOR) {
150  return true;
151  }
152  }
153  return false;
154 }
155 
156 /*!
157  * Return true if the current map-plane cell contains the grounding line, false otherwise.
158  *
159  * This is used to determine whether to use more quadrature points to estimate integrals
160  * over the bottom face of the basal element.
161  *
162  * The code takes advantage of the ordering of element nodes: lower 4 first, then upper 4.
163  * This means that we can loop over the first 4 nodes and ignore the other 4.
164  */
165 bool Blatter::grounding_line(const double *F) {
166 
167  // number of nodes per map-plane cell
168  int N = 4;
169 
170  bool
171  grounded = false,
172  floating = false;
173 
174  for (int n = 0; n < N; ++n) {
175  if (F[n] <= 0.0) {
176  grounded = true;
177  } else {
178  floating = true;
179  }
180  }
181 
182  return grounded and floating;
183 }
184 
185 /*!
186  * Return true if the current vertical face is partially submerged.
187  *
188  * This is used to determine whether to use more quadrature points to estimate integrals
189  * over this face when computing lateral boundary conditions.
190  */
191 bool Blatter::partially_submerged_face(int face, const double *z, const double *sea_level) {
192  const auto *nodes = fem::q13d::incident_nodes[face];
193 
194  // number of nodes per face
195  int N = 4;
196 
197  bool
198  above = false,
199  below = false;
200 
201  for (int n = 0; n < N; ++n) {
202  int k = nodes[n];
203  if (z[k] > sea_level[k]) {
204  above = true;
205  } else {
206  below = true;
207  }
208  }
209 
210  return above and below;
211 }
212 
213 /*!
214  * Return true if the current face is a part of the marine ice boundary (i.e. at a
215  * partially-submerged vertical cliff), false otherwise.
216  *
217  * A face is a part of the marine boundary if all four nodes are boundary nodes *and* at
218  * least one map-plane location has bottom elevation below sea level (floatation level).
219  *
220  * If a node is *both* a boundary and a Dirichlet node (this may happen), then we treat it
221  * as a boundary node here: element.add_contribution() will do the right thing in this
222  * case.
223  */
225  const int *node_type,
226  const double *ice_bottom,
227  const double *sea_level) {
228  const auto *nodes = fem::q13d::incident_nodes[face];
229 
230  // number of nodes per face
231  int N = 4;
232 
233  // exclude faces that contain at least one node that is not a part of the boundary
234  for (int n = 0; n < N; ++n) {
235  if (not (node_type[nodes[n]] == NODE_BOUNDARY)) {
236  return false;
237  }
238  }
239 
240  // This face is a part of the lateral boundary. Now we need to check if ice_bottom is
241  // below sea_level at one of the nodes of this face.
242  for (int n = 0; n < N; ++n) {
243  if (ice_bottom[nodes[n]] < sea_level[nodes[n]]) {
244  return true;
245  }
246  }
247  return false;
248 }
249 
250 /*!
251  * Allocate the Blatter-Pattyn stress balance solver
252  *
253  * @param[in] grid PISM's grid.
254  * @param[in] Mz number of vertical levels
255  * @param[in] n_levels maximum number of grid levels to use
256  * @param[in] coarsening_factor grid coarsening factor
257  */
258 Blatter::Blatter(IceGrid::ConstPtr grid, int Mz, int coarsening_factor)
259  : ShallowStressBalance(grid),
260  m_parameters(grid, "bp_input_parameters", WITH_GHOSTS),
261  m_face4(grid->dx(), grid->dy(), fem::Q1Quadrature4()), // 4-point Gaussian quadrature
262  m_face100(grid->dx(), grid->dy(), fem::Q1QuadratureN(10)) // 100-point quadrature for grounding lines
263 {
264 
265  assert(m_face4.n_pts() <= m_Nq);
266  assert(m_face100.n_pts() <= m_Nq);
267 
268  auto pism_da = grid->get_dm(1, 0);
269 
270  int ierr = setup(*pism_da, grid->periodicity(), Mz, coarsening_factor, "bp_");
271  if (ierr != 0) {
273  "Failed to allocate a Blatter solver instance");
274  }
275 
276  {
277  std::vector<double> sigma(Mz);
278  double dz = 1.0 / (Mz - 1.0);
279  for (int i = 0; i < Mz; ++i) {
280  sigma[i] = i * dz;
281  }
282  sigma.back() = 1.0;
283 
284  m_u_sigma.reset(new IceModelVec3(grid, "uvel_sigma", WITHOUT_GHOSTS, sigma));
285  m_u_sigma->set_attrs("diagnostic",
286  "u velocity component on the sigma grid",
287  "m s-1", "m s-1", "", 0);
288 
289  m_v_sigma.reset(new IceModelVec3(grid, "vvel_sigma", WITHOUT_GHOSTS, sigma));
290  m_v_sigma->set_attrs("diagnostic",
291  "v velocity component on the sigma grid",
292  "m s-1", "m s-1", "", 0);
293 
294  std::map<std::string,std::string> z_attrs =
295  {{"axis", "Z"},
296  {"long_name", "scaled Z-coordinate in the ice (z_base=0, z_surface=1)"},
297  {"units", "1"},
298  {"positive", "up"}};
299 
300  m_u_sigma->metadata(0).z().set_name("z_sigma");
301  m_v_sigma->metadata(0).z().set_name("z_sigma");
302 
303  for (const auto &z_attr : z_attrs) {
304  m_u_sigma->metadata(0).z().set_string(z_attr.first, z_attr.second);
305  m_v_sigma->metadata(0).z().set_string(z_attr.first, z_attr.second);
306  }
307 
308  }
309 
310  {
311  rheology::FlowLawFactory ice_factory("stress_balance.blatter.", m_config, m_EC);
312  ice_factory.remove(ICE_GOLDSBY_KOHLSTEDT);
313  m_flow_law = ice_factory.create();
314  }
315 
316  double g = m_config->get_number("constants.standard_gravity");
317 
318  m_rho_ice_g = m_config->get_number("constants.ice.density") * g;
319  m_rho_ocean_g = m_config->get_number("constants.sea_water.density") * g;
320 
321  m_eta_transform = m_config->get_flag("stress_balance.blatter.use_eta_transform");
322 
323  m_glen_exponent = m_flow_law->exponent();
324 
325  double E = m_config->get_number("stress_balance.blatter.enhancement_factor");
326  m_E_viscosity = std::pow(E, -1.0 / m_glen_exponent);
327 
328  double softness = m_config->get_number("flow_law.isothermal_Glen.ice_softness"),
329  hardness = std::pow(softness, -1.0 / m_glen_exponent);
330 
331  m_scaling = m_rho_ice_g * hardness;
332 }
333 
334 /*!
335  * Restrict 2D and 3D model parameters from a fine grid to a coarse grid.
336  *
337  * Re-compute node types from geometry.
338  *
339  * This hook is called every time SNES needs to update coarse-grid data.
340  *
341  * FIXME: parameters restricted by this hook do not change from one SNES iteration to the
342  * next, so we can return early after the first one.
343  */
344 static PetscErrorCode blatter_restriction_hook(DM fine,
345  Mat mrestrict, Vec rscale, Mat inject,
346  DM coarse, void *ctx) {
347  // Get rid of warnings about unused arguments
348  (void) mrestrict;
349  (void) rscale;
350  (void) inject;
351  (void) ctx;
352 
353  PetscErrorCode ierr = restrict_data(fine, coarse, "3D_DM"); CHKERRQ(ierr);
354 
355  return 0;
356 }
357 
358 static PetscErrorCode blatter_coarsening_hook(DM dm_fine, DM dm_coarse, void *ctx) {
359  PetscErrorCode ierr;
360 
361  int mg_levels = 1;
362  {
363  const char *prefix;
364  ierr = DMGetOptionsPrefix(dm_fine, &prefix); CHKERRQ(ierr);
365  auto option = pism::printf("-%spc_mg_levels", prefix);
366  mg_levels = options::Integer(option, "", 1);
367  }
368 
369  ierr = setup_level(dm_coarse, mg_levels); CHKERRQ(ierr);
370 
371  ierr = DMCoarsenHookAdd(dm_coarse, blatter_coarsening_hook, blatter_restriction_hook, ctx); CHKERRQ(ierr);
372 
373  // 3D
374  ierr = create_restriction(dm_fine, dm_coarse, "3D_DM"); CHKERRQ(ierr);
375 
376  return 0;
377 }
378 
379 /*!
380  * Allocates the 3D DM, the corresponding solution vector, and the SNES solver.
381  */
382 PetscErrorCode Blatter::setup(DM pism_da, Periodicity periodicity, int Mz,
383  int coarsening_factor,
384  const std::string &prefix) {
385  MPI_Comm comm;
386  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pism_da, &comm); CHKERRQ(ierr);
387 
388  // FIXME: add the ability to add a prefix to the option prefix. We need this to be able
389  // to run more than one instance of PISM in parallel.
390  auto option = pism::printf("-%spc_mg_levels", prefix.c_str());
391  int mg_levels = options::Integer(option, "", 1);
392 
393  // Check compatibility of Mz, mg_levels, and the coarsening_factor and stop if they are
394  // not compatible.
395  //
396  // We assume that the user also set "-bp_pc_type mg".
397  {
398  int c = coarsening_factor;
399  int M = mg_levels;
400  int mz = Mz;
401  while (M > 1) {
402  // Note: integer division
403  if (((mz - 1) / c) * c != mz - 1) {
404  int N = std::pow(c, (int)mg_levels - 1);
405  auto message = pism::printf("Blatter stress balance solver: settings\n"
406  "stress_balance.blatter.Mz = %d,\n"
407  "stress_balance.blatter.coarsening_factor = %d,\n"
408  "and '%s %d' are not compatible.\n"
409  "To use N = %d multigrid levels with the coarsening factor C = %d\n"
410  "stress_balance.blatter.Mz has to be equal to A * C^(N - 1) + 1\n"
411  "for some positive integer A, e.g. %d, %d, %d, ...",
412  Mz, c, option.c_str(), mg_levels, mg_levels, c,
413  N + 1, 2*N + 1, 3*N + 1);
414  throw RuntimeError(PISM_ERROR_LOCATION, message);
415  }
416  mz = (mz - 1) / c + 1;
417  M -= 1;
418  }
419  }
420 
421  // DM
422  //
423  // Note: in the PISM's DA pism_da PETSc's and PISM's meaning of x and y are the same.
424  {
425  PetscInt Mx, My;
426  PetscInt mx, my;
427  PetscInt dims;
428 
429  ierr = DMDAGetInfo(pism_da,
430  &dims, // dimensions
431  &Mx, &My, NULL, // grid size
432  &mx, &my, NULL, // numbers of processors in each direction
433  NULL, // number of degrees of freedom
434  NULL, // stencil width
435  NULL, NULL, NULL, // types of ghost nodes at the boundary
436  NULL); // stencil type
437  CHKERRQ(ierr);
438 
439  assert(dims == 2);
440 
441  const PetscInt
442  *lx = NULL,
443  *ly = NULL;
444  ierr = DMDAGetOwnershipRanges(pism_da, &lx, &ly, NULL); CHKERRQ(ierr);
445 
446  PetscInt
447  mz = 1,
448  dof = 2,
449  stencil_width = 1;
450 
451  DMBoundaryType
452  bx = (periodicity & X_PERIODIC) != 0 ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
453  by = (periodicity & Y_PERIODIC) != 0 ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
454  bz = DM_BOUNDARY_NONE;
455 
456  ierr = DMDACreate3d(comm,
457  bz, bx, by, // STORAGE_ORDER
458  DMDA_STENCIL_BOX,
459  Mz, Mx, My, // STORAGE_ORDER
460  mz, mx, my, // STORAGE_ORDER
461  dof,
462  stencil_width,
463  NULL, lx, ly, // STORAGE_ORDER
464  m_da.rawptr()); CHKERRQ(ierr);
465 
466  ierr = DMSetOptionsPrefix(m_da, prefix.c_str()); CHKERRQ(ierr);
467 
468  // semi-coarsening: coarsen in the vertical direction only
469  ierr = DMDASetRefinementFactor(m_da, coarsening_factor, 1, 1); CHKERRQ(ierr); // STORAGE_ORDER
470 
471  ierr = DMSetFromOptions(m_da); CHKERRQ(ierr);
472 
473  ierr = DMSetUp(m_da); CHKERRQ(ierr);
474 
475  // set up 3D parameter storage
476  ierr = setup_level(m_da, mg_levels); CHKERRQ(ierr);
477 
478  // tell PETSc how to coarsen this grid and how to restrict data to a coarser grid
479  ierr = DMCoarsenHookAdd(m_da, blatter_coarsening_hook, blatter_restriction_hook, NULL);
480  CHKERRQ(ierr);
481  }
482 
483  // Vec
484  {
485  ierr = DMCreateGlobalVector(m_da, m_x.rawptr()); CHKERRQ(ierr);
486 
487  ierr = VecSetOptionsPrefix(m_x, prefix.c_str()); CHKERRQ(ierr);
488 
489  ierr = VecSetFromOptions(m_x); CHKERRQ(ierr);
490 
491  ierr = VecDuplicate(m_x, m_x_old.rawptr()); CHKERRQ(ierr);
492  }
493 
494  // SNES
495  {
496  ierr = SNESCreate(comm, m_snes.rawptr()); CHKERRQ(ierr);
497 
498  ierr = SNESSetOptionsPrefix(m_snes, prefix.c_str()); CHKERRQ(ierr);
499 
500  ierr = SNESSetDM(m_snes, m_da); CHKERRQ(ierr);
501 
502  ierr = DMDASNESSetFunctionLocal(m_da, INSERT_VALUES,
503  (DMDASNESFunction)function_callback,
504  this); CHKERRQ(ierr);
505 
506  ierr = DMDASNESSetJacobianLocal(m_da,
507  (DMDASNESJacobian)jacobian_callback,
508  this); CHKERRQ(ierr);
509 
510  ierr = SNESSetFromOptions(m_snes); CHKERRQ(ierr);
511 
512 
513  PetscBool ksp_use_ew = PETSC_FALSE;
514  ierr = SNESKSPGetUseEW(m_snes, &ksp_use_ew); CHKERRQ(ierr);
515  m_ksp_use_ew = (ksp_use_ew != 0U);
516  }
517 
518  return 0;
519 }
520 
521 /*!
522  * Set 2D parameters on the finest grid.
523  */
524 void Blatter::init_2d_parameters(const Inputs &inputs) {
525 
526  double
527  ice_density = m_config->get_number("constants.ice.density"),
528  water_density = m_config->get_number("constants.sea_water.density"),
529  alpha = ice_density / water_density;
530 
531  const IceModelVec2S
532  &tauc = *inputs.basal_yield_stress,
533  &H = inputs.geometry->ice_thickness,
534  &b = inputs.geometry->bed_elevation,
535  &sea_level = inputs.geometry->sea_level_elevation;
536 
537  {
538  IceModelVec::AccessList list{&tauc, &H, &b, &sea_level, &m_parameters};
539 
540  for (Points p(*m_grid); p; p.next()) {
541  const int i = p.i(), j = p.j();
542 
543  double
544  b_grounded = b(i, j),
545  b_floating = sea_level(i, j) - alpha * H(i, j),
546  s_grounded = b(i, j) + H(i, j),
547  s_floating = sea_level(i, j) + (1.0 - alpha) * H(i, j);
548 
549  m_parameters(i, j).tauc = tauc(i, j);
550  m_parameters(i, j).thickness = H(i, j);
551  m_parameters(i, j).sea_level = sea_level(i, j);
552  m_parameters(i, j).bed = std::max(b_grounded, b_floating);
553  m_parameters(i, j).node_type = NODE_EXTERIOR;
554  m_parameters(i, j).floatation = s_floating - s_grounded;
555  }
556  }
557 
558  // update ghosts here: the call to compute_node_type() uses ghosts of ice thickness
559  m_parameters.update_ghosts();
560 
561  double H_min = m_config->get_number("stress_balance.ice_free_thickness_standard");
562  compute_node_type(H_min);
563 
564  // update ghosts of node types stored in m_parameters
565  m_parameters.update_ghosts();
566 }
567 
568 /*!
569  * Set 3D parameters on the finest grid.
570  */
571 void Blatter::init_ice_hardness(const Inputs &inputs, const petsc::DM &da) {
572 
573  const auto *enthalpy = inputs.enthalpy;
574  // PISM's vertical grid:
575  const auto &zlevels = enthalpy->levels();
576  auto Mz = zlevels.size();
577 
578  // solver's vertical grid:
579  int Mz_sigma = 0;
580  {
581  DMDALocalInfo info;
582  int ierr = DMDAGetLocalInfo(da, &info); PISM_CHK(ierr, "DMDAGetLocalInfo");
583  info = grid_transpose(info);
584  Mz_sigma = info.mz;
585  }
586 
587  const auto &ice_thickness = inputs.geometry->ice_thickness;
588  DataAccess<double***> hardness(da, 3, NOT_GHOSTED);
589 
590  IceModelVec::AccessList list{enthalpy, &ice_thickness};
591 
592  for (Points p(*m_grid); p; p.next()) {
593  const int i = p.i(), j = p.j();
594 
595  double H = ice_thickness(i, j);
596 
597  const double *E = enthalpy->get_column(i, j);
598 
599  for (int k = 0; k < Mz_sigma; ++k) {
600  double
601  z = grid_z(0.0, H, Mz_sigma, k),
602  depth = H - z,
603  pressure = m_EC->pressure(depth),
604  E_local = 0.0;
605 
606  auto k0 = m_grid->kBelowHeight(z);
607 
608  if (k0 + 1 < Mz) {
609  double lambda = (z - zlevels[k0]) / (zlevels[k0 + 1] - zlevels[k0]);
610 
611  E_local = (1.0 - lambda) * E[k0] + lambda * E[k0 + 1];
612  } else {
613  E_local = E[Mz - 1];
614  }
615 
616  hardness[j][i][k] = m_flow_law->hardness(E_local, pressure); // STORAGE_ORDER
617 
618  } // end of the loop over sigma levels
619  } // end of the loop over grid points
620 }
621 
622 /*!
623  * Get values of 2D parameters at element nodes.
624  *
625  * This method is re-implemented by derived classes that use periodic boundary conditions.
626  */
628  Parameters **P,
629  int i,
630  int j,
631  int *node_type,
632  double *bottom_elevation,
633  double *ice_thickness,
634  double *surface_elevation,
635  double *sea_level) const {
636 
637  for (int n = 0; n < fem::q13d::n_chi; ++n) {
638  auto I = element.local_to_global(i, j, 0, n);
639 
640  auto p = P[I.j][I.i];
641 
642  node_type[n] = static_cast<int>(p.node_type);
643  bottom_elevation[n] = p.bed;
644  ice_thickness[n] = p.thickness;
645 
646  if (surface_elevation != nullptr) {
647  surface_elevation[n] = p.bed + p.thickness;
648  }
649 
650  if (sea_level != nullptr) {
651  sea_level[n] = p.sea_level;
652  }
653  }
654 }
655 
657  m_log->message(2, "* Initializing the Blatter stress balance...\n");
658 
660 
661  if (opts.type == INIT_RESTART) {
662  File input_file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
663  bool u_sigma_found = input_file.find_variable("uvel_sigma");
664  bool v_sigma_found = input_file.find_variable("vvel_sigma");
665  unsigned int start = input_file.nrecords() - 1;
666 
667  if (u_sigma_found and v_sigma_found) {
668  m_log->message(3, "Reading uvel_sigma and vvel_sigma...\n");
669 
670  m_u_sigma->read(input_file, start);
671  m_v_sigma->read(input_file, start);
672 
674  } else {
676  "uvel_sigma and vvel_sigma not found");
677  }
678  } else {
679  int ierr = VecSet(m_x, 0.0); PISM_CHK(ierr, "VecSet");
680  }
681 }
682 
683 void Blatter::define_model_state_impl(const File &output) const {
684  m_u_sigma->define(output);
685  m_v_sigma->define(output);
686 }
687 
688 void Blatter::write_model_state_impl(const File &output) const {
689  m_u_sigma->write(output);
690  m_v_sigma->write(output);
691 }
692 
694 
695  DMDALocalInfo info;
696  int ierr = DMDAGetLocalInfo(m_da, &info); PISM_CHK(ierr, "DMDAGetLocalInfo");
697  info = grid_transpose(info);
698 
699  fem::Q1Element2 E(info, 1.0, 1.0, fem::Q1Quadrature1());
700 
702 
703  double R_min = 1e16, R_max = 0.0, R_avg = 0.0;
704  double dxy = std::max(m_grid->dx(), m_grid->dy());
705  double n_cells = 0.0;
706 
708  for (int j = info.ys; j < info.ys + info.ym - 1; j++) {
709  for (int i = info.xs; i < info.xs + info.xm - 1; i++) {
710 
711  E.reset(i, j);
712 
713  E.nodal_values(m_parameters.array(), P);
714 
715  int node_type[4];
716  for (int k = 0; k < 4; ++k) {
717  node_type[k] = static_cast<int>(P[k].node_type);
718  }
719 
720  if (exterior_element(node_type)) {
721  continue;
722  }
723 
724  n_cells += 1.0;
725 
726  double dz_max = 0.0;
727  for (int k = 0; k < 4; ++k) {
728  dz_max = std::max(P[k].thickness / (info.mz - 1), dz_max);
729  }
730  double R = dz_max / dxy;
731 
732  R_min = std::min(R, R_min);
733  R_max = std::max(R, R_max);
734  R_avg += R;
735  }
736  }
737 
738  n_cells = GlobalSum(m_grid->com, n_cells);
739  R_avg = GlobalSum(m_grid->com, R_avg);
740  R_avg /= std::max(n_cells, 1.0);
741 
742  R_min = GlobalMin(m_grid->com, R_min);
743  R_max = GlobalMax(m_grid->com, R_max);
744 
745  m_log->message(2,
746  "Blatter solver: %d * (%d - 1) = %d active elements\n",
747  (int)n_cells, (int)info.mz, (int)(n_cells * (info.mz - 1)));
748 
749  if (n_cells > 0) {
750  m_log->message(2,
751  " Vertical spacing (m): min = %3.3f, max = %3.3f, avg = %3.3f\n",
752  R_min * dxy, R_max * dxy, R_avg * dxy);
753  m_log->message(2,
754  " Aspect ratios: min = %3.3f, max = %3.3f, avg = %3.3f, max/min = %3.3f\n",
755  R_min, R_max, R_avg, R_max / R_min);
756  }
757 }
758 
759 /*!
760  * Runs the solver and extracts iteration counts.
761  */
763  PetscErrorCode ierr;
764  SolutionInfo result;
765 
766  // Solve the system:
767  ierr = SNESSolve(m_snes, NULL, m_x); PISM_CHK(ierr, "SNESSolve");
768 
769  ierr = SNESGetConvergedReason(m_snes, &result.snes_reason);
770  PISM_CHK(ierr, "SNESGetConvergedReason");
771 
772  ierr = SNESGetIterationNumber(m_snes, &result.snes_it);
773  PISM_CHK(ierr, "SNESGetIterationNumber");
774 
775  ierr = SNESGetLinearSolveIterations(m_snes, &result.ksp_it);
776  PISM_CHK(ierr, "SNESGetLinearSolveIterations");
777 
778  KSP ksp;
779  ierr = SNESGetKSP(m_snes, &ksp);
780  PISM_CHK(ierr, "SNESGetKSP");
781 
782  ierr = KSPGetConvergedReason(ksp, &result.ksp_reason);
783  PISM_CHK(ierr, "KSPGetConvergedReason");
784 
785  PC pc;
786  ierr = KSPGetPC(ksp, &pc);
787  PISM_CHK(ierr, "KSPGetPC");
788 
789  PCType pc_type;
790  ierr = PCGetType(pc, &pc_type);
791  PISM_CHK(ierr, "PCGetType");
792 
793  if (std::string(pc_type) == PCMG) {
794  KSP coarse_ksp;
795  ierr = PCMGGetCoarseSolve(pc, &coarse_ksp);
796  PISM_CHK(ierr, "PCMGGetCoarseSolve");
797 
798  ierr = KSPGetIterationNumber(coarse_ksp, &result.mg_coarse_ksp_it);
799  PISM_CHK(ierr, "KSPGetIterationNumber");
800  } else {
801  result.mg_coarse_ksp_it = 0;
802  }
803 
804  return result;
805 }
806 
808  PetscErrorCode ierr;
809 
810  SolutionInfo info;
811  // total number of SNES and KSP iterations
812  int snes_total_it = 0;
813  int ksp_total_it = 0;
814 
815  // maximum number of continuation steps (input)
816  int Nc = 50;
817 
818  double
819  schoofLen = m_config->get_number("flow_law.Schoof_regularizing_length", "m"),
820  schoofVel = m_config->get_number("flow_law.Schoof_regularizing_velocity", "m second-1"),
821  // desired regularization parameter
822  eps = pow(schoofVel / schoofLen, 2.0),
823  // gamma is a number such that 10^gamma <= eps. It is used to convert lambda in [0, 1] to eps_n
824  gamma = std::floor(std::log10(eps)),
825  // starting value of lambda (input)
826  lambda_min = 0.0,
827  // Final value of lambda (fixed)
828  lambda_max = 1.0,
829  // Minimum step length (input)
830  delta_min = 0.01,
831  // Maximum step length (input)
832  delta_max = 0.2,
833  // Initial increment of lambda (input)
834  delta0 = 0.05,
835  // "Aggressiveness" of the step increase, a non-negative number (input). The effect of
836  // this parameter is linked to the choice of -bp_snes_max_it obtained below.
837  A = 0.25;
838 
839  PetscInt snes_max_it = 0;
840  ierr = SNESGetTolerances(m_snes, NULL, NULL, NULL, &snes_max_it, NULL);
841  PISM_CHK(ierr, "SNESGetTolerances");
842 
843  double lambda = lambda_min;
844  double delta = delta0;
845 
846  // Use the zero initial guess to start
847  ierr = VecSet(m_x, 0.0); PISM_CHK(ierr, "VecSet");
848  ierr = VecSet(m_x_old, 0.0); PISM_CHK(ierr, "VecSet");
849 
850  m_log->message(2,
851  "Blatter solver: Starting parameter continuation with lambda = %f\n",
852  lambda);
853 
854  for (int N = 0; N < Nc; ++N) {
855  // Set the regularization parameter:
856  m_viscosity_eps = std::max(std::pow(10.0, lambda * gamma), eps);
857 
858  m_log->message(2, "Blatter solver: step %d with lambda = %f, eps = %e\n",
859  N, lambda, m_viscosity_eps);
860 
861  // Solve the system:
862 
863  info = solve();
864  snes_total_it += info.snes_it;
865  ksp_total_it += info.ksp_it;
866 
867  // report number of iterations for this continuation step
868  m_log->message(2, "Blatter solver continuation step #%d: %s\n"
869  " lambda = %f, eps = %e\n"
870  " SNES: %d, KSP: %d\n",
871  N, SNESConvergedReasons[info.snes_reason], lambda, m_viscosity_eps,
872  (int)info.snes_it, (int)info.ksp_it);
873  if (info.mg_coarse_ksp_it > 0) {
874  m_log->message(2,
875  " Coarse MG level KSP (last iteration): %d\n",
876  (int)info.mg_coarse_ksp_it);
877  }
878 
879  if (info.snes_reason > 0) {
880  // converged
881 
882  if (m_viscosity_eps <= eps) {
883  // ... while solving the desired (not overregularized) problem
884  info.snes_it = snes_total_it;
885  info.ksp_it = ksp_total_it;
886  return info;
887  }
888 
889  // Store the solution as the "old" initial guess we may need to revert to
890  ierr = VecCopy(m_x, m_x_old); PISM_CHK(ierr, "VecCopy");
891 
892  if (N > 0) {
893  // Adjust delta using the formula from LOCA (equation 2.8 in Salinger2002
894  // corrected using the code in Trilinos).
895  double F = (snes_max_it - info.snes_it) / (double)snes_max_it;
896  delta *= 1.0 + A * F * F;
897  }
898 
899  delta = std::min(delta, delta_max);
900 
901  // Ensure that delta does not take us past lambda_max
902  if (lambda + delta > lambda_max) {
903  delta = lambda_max - lambda;
904  }
905 
906  m_log->message(2, " Advancing lambda from %f to %f (delta = %f)\n",
907  lambda, lambda + delta, delta);
908 
909  lambda += delta;
910  } else if (info.snes_reason == SNES_DIVERGED_LINE_SEARCH or
911  info.snes_reason == SNES_DIVERGED_MAX_IT) {
912  // a continuation step failed
913 
914  // restore the previous initial guess
915  ierr = VecCopy(m_x_old, m_x); PISM_CHK(ierr, "VecCopy");
916 
917  // revert lambda to the previous value
918  lambda -= delta;
919 
920  if (lambda < lambda_min) {
921  m_log->message(2, "Blatter solver: Parameter continuation failed at step 1\n");
922  return info;
923  }
924 
925  if (std::fabs(delta - delta_min) < 1e-6) {
926  m_log->message(2, "Blatter solver: cannot reduce the continuation step\n");
927  return info;
928  }
929 
930  // reduce the step size
931  delta *= 0.5;
932 
933  delta = pism::clip(delta, delta_min, delta_max);
934 
935  lambda += delta;
936  // Note that this delta will not take us past lambda_max because the original
937  // delta satisfies lambda + delta <= lambda_max.
938 
939  m_log->message(2, " Back-tracking to lambda = %f using delta = %f\n",
940  lambda, delta);
941  } else {
942  return info;
943  }
944  }
945 
946  m_log->message(2, "Blatter solver failed after %d parameter continuation steps\n", Nc);
947 
948  return info;
949 }
950 
951 /*!
952  * Disable the Eisenstat-Walker method of setting KSP tolerances.
953  */
954 static void disable_ew(::SNES snes, double rtol) {
955  PetscErrorCode ierr;
956 
957  ierr = SNESKSPSetUseEW(snes, PETSC_FALSE);
958  PISM_CHK(ierr, "SNESKSPSetUseEW");
959 
960  KSP ksp;
961  ierr = SNESGetKSP(snes, &ksp);
962  PISM_CHK(ierr, "SNESGetKSP");
963 
964  ierr = KSPSetTolerances(ksp, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
965  PISM_CHK(ierr, "KSPSetTolerances");
966 }
967 
968 static void enable_ew(::SNES snes) {
969  PetscErrorCode ierr;
970  // restore the old EW setting
971  ierr = SNESKSPSetUseEW(snes, PETSC_TRUE); PISM_CHK(ierr, "SNESKSPSetUseEW");
972 }
973 
974 void Blatter::update(const Inputs &inputs, bool full_update) {
975  PetscErrorCode ierr;
976  (void) full_update;
977 
978  {
979  double
980  schoofLen = m_config->get_number("flow_law.Schoof_regularizing_length", "m"),
981  schoofVel = m_config->get_number("flow_law.Schoof_regularizing_velocity", "m second-1");
982 
983  m_viscosity_eps = pow(schoofVel / schoofLen, 2.0);
984 
985  }
986 
987  init_2d_parameters(inputs);
988  init_ice_hardness(inputs, m_da);
989 
991 
992  // Store the "old" initial guess: it may be needed to re-try.
993  ierr = VecCopy(m_x, m_x_old); PISM_CHK(ierr, "VecCopy");
994 
995  SolutionInfo info;
996  int snes_total_it = 0;
997  int ksp_total_it = 0;
998 
999  double ksp_rtol = 1e-5;
1000 
1001  double norm = 0.0;
1002  {
1003  ierr = VecNorm(m_x, NORM_INFINITY, &norm); PISM_CHK(ierr, "VecNorm");
1004  }
1005 
1006  // First attempt
1007  {
1008  if (m_ksp_use_ew and norm == 0.0) {
1009  m_log->message(2,
1010  "Blatter solver: zero initial guess\n"
1011  " Disabling the Eisenstat-Walker method of adjusting solver tolerances\n");
1012 
1013  disable_ew(m_snes, ksp_rtol);
1014  info = solve();
1015  enable_ew(m_snes);
1016  } else {
1017  info = solve();
1018  }
1019  snes_total_it += info.snes_it;
1020  ksp_total_it += info.ksp_it;
1021 
1022  if (info.snes_reason > 0) {
1023  goto bp_done;
1024  }
1025  m_log->message(2, "Blatter solver: %s\n", SNESConvergedReasons[info.snes_reason]);
1026  }
1027 
1028  if (m_ksp_use_ew and norm > 0.0) {
1029  m_log->message(2," Trying without the Eisenstat-Walker method of adjusting solver tolerances\n");
1030 
1031  // restore the previous initial guess
1032  if (not (info.snes_reason == SNES_DIVERGED_LINE_SEARCH or
1033  info.snes_reason == SNES_DIVERGED_MAX_IT)) {
1034  ierr = VecCopy(m_x_old, m_x); PISM_CHK(ierr, "VecCopy");
1035  } else {
1036  // We *keep* the current values in m_x if the line search failed after a few
1037  // iterations or if the solver took too many iterations: no need to discard the
1038  // progress it made before failing.
1039  }
1040 
1041  {
1042  disable_ew(m_snes, ksp_rtol);
1043  info = solve();
1044  enable_ew(m_snes);
1045  }
1046 
1047  snes_total_it += info.snes_it;
1048  ksp_total_it += info.ksp_it;
1049 
1050  if (info.snes_reason > 0) {
1051  goto bp_done;
1052  }
1053  m_log->message(2, "Blatter solver: %s\n", SNESConvergedReasons[info.snes_reason]);
1054  }
1055 
1056  // try using parameter continuation
1057  {
1058  info = parameter_continuation();
1059  snes_total_it += info.snes_it;
1060  ksp_total_it += info.ksp_it;
1061  if (info.snes_reason > 0) {
1062  goto bp_done;
1063  }
1064  m_log->message(2, "Blatter solver: %s\n", SNESConvergedReasons[info.snes_reason]);
1065  }
1066 
1067  throw RuntimeError(PISM_ERROR_LOCATION, "Blatter solver failed");
1068 
1069  bp_done:
1070  // report the total number of iterations
1071  m_log->message(2,
1072  "Blatter solver: %s. Done.\n"
1073  " SNES: %d, KSP: %d\n",
1074  SNESConvergedReasons[info.snes_reason],
1075  snes_total_it, ksp_total_it);
1076  if (info.mg_coarse_ksp_it > 0) {
1077  m_log->message(2,
1078  " Level 0 KSP (last iteration): %d\n",
1079  (int)info.mg_coarse_ksp_it);
1080  }
1081 
1082  // put basal velocity in m_velocity to use it in the next call
1084 
1086  inputs.geometry->cell_type,
1088 
1090 
1091  // copy the solution from m_x to m_u_sigma, m_v_sigma for re-starting
1092  copy_solution();
1093 }
1094 
1096  Vector2 ***x = nullptr;
1097  int ierr = DMDAVecGetArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecGetArray");
1098 
1099  int Mz = m_u_sigma->levels().size();
1100 
1101  IceModelVec::AccessList list{m_u_sigma.get(), m_v_sigma.get()};
1102 
1103  for (Points p(*m_grid); p; p.next()) {
1104  const int i = p.i(), j = p.j();
1105 
1106  auto *u = m_u_sigma->get_column(i, j);
1107  auto *v = m_v_sigma->get_column(i, j);
1108 
1109  for (int k = 0; k < Mz; ++k) {
1110  u[k] = x[j][i][k].u; // STORAGE_ORDER
1111  v[k] = x[j][i][k].v; // STORAGE_ORDER
1112  }
1113  }
1114 
1115  ierr = DMDAVecRestoreArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecRestoreArray");
1116 }
1117 
1119  Vector2 ***x = nullptr;
1120  int ierr = DMDAVecGetArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecGetArray");
1121 
1122  IceModelVec::AccessList list{&result};
1123 
1124  for (Points p(*m_grid); p; p.next()) {
1125  const int i = p.i(), j = p.j();
1126 
1127  result(i, j) = x[j][i][0]; // STORAGE_ORDER
1128  }
1129 
1130  ierr = DMDAVecRestoreArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecRestoreArray");
1131 }
1132 
1133 
1135  const IceModelVec3 &v_sigma) {
1136  Vector2 ***x = nullptr;
1137  int ierr = DMDAVecGetArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecGetArray");
1138 
1139  int Mz = m_u_sigma->levels().size();
1140 
1141  IceModelVec::AccessList list{&u_sigma, &v_sigma};
1142 
1143  for (Points p(*m_grid); p; p.next()) {
1144  const int i = p.i(), j = p.j();
1145 
1146  const auto *u = u_sigma.get_column(i, j);
1147  const auto *v = v_sigma.get_column(i, j);
1148 
1149  for (int k = 0; k < Mz; ++k) {
1150  x[j][i][k].u = u[k]; // STORAGE_ORDER
1151  x[j][i][k].v = v[k]; // STORAGE_ORDER
1152  }
1153  }
1154 
1155  ierr = DMDAVecRestoreArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecRestoreArray");
1156 }
1157 
1159  PetscErrorCode ierr;
1160 
1161  Vector2 ***x = nullptr;
1162  ierr = DMDAVecGetArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecGetArray");
1163 
1164  int Mz = m_u_sigma->levels().size();
1165 
1166  IceModelVec::AccessList list{&result, &m_parameters};
1167 
1168  for (Points p(*m_grid); p; p.next()) {
1169  const int i = p.i(), j = p.j();
1170 
1171  double H = m_parameters(i, j).thickness;
1172 
1173  Vector2 V(0.0, 0.0);
1174 
1175  if (H > 0.0) {
1176  // use trapezoid rule to compute the column average
1177  double dz = H / (Mz - 1);
1178  for (int k = 0; k < Mz - 1; ++k) {
1179  V += x[j][i][k] + x[j][i][k + 1]; // STORAGE_ORDER
1180  }
1181  V *= (0.5 * dz) / H;
1182  }
1183 
1184  result(i, j) = V;
1185  }
1186 
1187  ierr = DMDAVecRestoreArray(m_da, m_x, &x); PISM_CHK(ierr, "DMDAVecRestoreArray");
1188 
1189  result.update_ghosts();
1190 }
1191 
1192 
1194  return m_u_sigma;
1195 }
1196 
1198  return m_v_sigma;
1199 }
1200 
1201 } // end of namespace stressbalance
1202 } // end of namespace pism
#define ICE_GOLDSBY_KOHLSTEDT
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
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
Definition: File.cc:364
unsigned int nrecords() const
Get the number of records. Uses the length of an unlimited dimension.
Definition: File.cc:316
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S sea_level_elevation
Definition: Geometry.hh:50
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
std::shared_ptr< IceModelVec3 > Ptr
Definition: iceModelVec.hh:423
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
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
std::vector< double > levels() const
Definition: iceModelVec.cc:136
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
void local_to_global(unsigned int k, int &i, int &j) const
Convert a local degree of freedom index k to a global degree of freedom index (i,j).
Definition: Element.hh:171
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
Definition: Element.cc:187
void nodal_values(const IceModelVec2Int &x_global, int *result) const
Get nodal values of an integer mask.
Definition: Element.cc:176
GlobalIndex local_to_global(int i, int j, int k, unsigned int n) const
Definition: Element.hh:298
Q1 element with sides parallel to X and Y axes.
Definition: Element.hh:257
int n_pts() const
Number of quadrature points.
Definition: Element.hh:404
3D Q1 element
Definition: Element.hh:351
The 1-point Gaussian quadrature on the square [-1,1]*[-1,1].
Definition: Quadrature.hh:89
std::shared_ptr< FlowLaw > create()
void remove(const std::string &name)
void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Blatter.cc:683
void compute_node_type(double min_thickness)
Definition: Blatter.cc:55
void set_initial_guess(const IceModelVec3 &u_sigma, const IceModelVec3 &v_sigma)
Definition: Blatter.cc:1134
void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Blatter.cc:688
virtual bool marine_boundary(int face, const int *node_type, const double *ice_bottom, const double *sea_level)
Definition: Blatter.cc:224
fem::Q1Element3Face m_face4
Definition: Blatter.hh:112
static bool partially_submerged_face(int face, const double *z, const double *sea_level)
Definition: Blatter.cc:191
IceModelVec3::Ptr m_u_sigma
Definition: Blatter.hh:68
PetscErrorCode setup(DM pism_da, Periodicity p, int Mz, int coarsening_factor, const std::string &prefix)
Definition: Blatter.cc:382
virtual void init_2d_parameters(const Inputs &inputs)
Definition: Blatter.cc:524
static PetscErrorCode function_callback(DMDALocalInfo *info, const Vector2 ***x, Vector2 ***f, Blatter *solver)
Definition: residual.cc:457
IceModelVec3::Ptr velocity_v_sigma() const
Definition: Blatter.cc:1197
IceModelVec3::Ptr velocity_u_sigma() const
Definition: Blatter.cc:1193
SolutionInfo parameter_continuation()
Definition: Blatter.cc:807
void update(const Inputs &inputs, bool)
Definition: Blatter.cc:974
IceModelVec3::Ptr m_v_sigma
Definition: Blatter.hh:68
void init_ice_hardness(const Inputs &inputs, const petsc::DM &da)
Definition: Blatter.cc:571
static const int m_Nq
Definition: Blatter.hh:105
fem::Q1Element3Face m_face100
Definition: Blatter.hh:113
void get_basal_velocity(IceModelVec2V &result)
Definition: Blatter.cc:1118
virtual Vector2 u_bc(double x, double y, double z) const
Definition: Blatter.cc:131
void compute_averaged_velocity(IceModelVec2V &result)
Definition: Blatter.cc:1158
static bool exterior_element(const int *node_type)
Definition: Blatter.cc:145
Blatter(IceGrid::ConstPtr grid, int Mz, int coarsening_factor)
Definition: Blatter.cc:258
static bool grounding_line(const double *F)
Definition: Blatter.cc:165
virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex &I)
Definition: Blatter.cc:123
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, const Vector2 ***x, Mat A, Mat J, Blatter *solver)
Definition: jacobian.cc:354
IceModelVec2< Parameters > m_parameters
Definition: Blatter.hh:79
virtual void nodal_parameter_values(const fem::Q1Element3 &element, Parameters **P, int i, int j, int *node_type, double *bottom, double *thickness, double *surface, double *sea_level) const
Definition: Blatter.cc:627
const IceModelVec3 * enthalpy
const IceModelVec2S * basal_yield_stress
std::shared_ptr< rheology::FlowLaw > m_flow_law
void compute_basal_frictional_heating(const IceModelVec2V &velocity, const IceModelVec2S &tauc, const IceModelVec2CellType &mask, IceModelVec2S &result) const
Compute the basal frictional heating.
Shallow stress balance (such as the SSA).
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
const int n_chi
Number of shape functions on a Q1 element.
Definition: FEM.hh:213
const unsigned int incident_nodes[n_faces][4]
Definition: FEM.hh:223
const int n_chi
Definition: FEM.hh:191
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:43
static void disable_ew(::SNES snes, double rtol)
Definition: Blatter.cc:954
static void enable_ew(::SNES snes)
Definition: Blatter.cc:968
static PetscErrorCode blatter_restriction_hook(DM fine, Mat mrestrict, Vec rscale, Mat inject, DM coarse, void *ctx)
Definition: Blatter.cc:344
static PetscErrorCode blatter_coarsening_hook(DM dm_fine, DM dm_coarse, void *ctx)
Definition: Blatter.cc:358
double grid_z(double b, double H, int Mz, int k)
PetscErrorCode setup_level(DM dm, int mg_levels)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static double F(double SL, double B, double H, double alpha)
@ NODE_BOUNDARY
Definition: node_types.hh:35
@ NODE_EXTERIOR
Definition: node_types.hh:36
@ NODE_INTERIOR
Definition: node_types.hh:34
@ NOT_GHOSTED
Definition: DataAccess.hh:30
Periodicity
Definition: IceGrid.hh:48
@ Y_PERIODIC
Definition: IceGrid.hh:48
@ X_PERIODIC
Definition: IceGrid.hh:48
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
Definition: Component.cc:45
static const double g
Definition: exactTestP.cc:39
DMDALocalInfo grid_transpose(const DMDALocalInfo &input)
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
@ INIT_RESTART
Definition: Component.hh:39
T clip(T x, T a, T b)
std::string printf(const char *format,...)
static const double k
Definition: exactTestP.cc:45
@ PISM_GUESS
Definition: IO_Flags.hh:41
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
PetscErrorCode restrict_data(DM fine, DM coarse, const char *dm_name)
Restrict model parameters from the fine grid onto the coarse grid.
void GlobalMin(MPI_Comm comm, double *local, double *result, int count)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
PetscErrorCode create_restriction(DM fine, DM coarse, const char *dm_name)
Create the restriction matrix.
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
const int I[]
Definition: ssafd_code.cc:24
InitializationType type
initialization type
Definition: Component.hh:44
std::string filename
name of the input file (if applicable)
Definition: Component.hh:46