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