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
SSAFEM.cc
Go to the documentation of this file.
1// Copyright (C) 2009--2025 Jed Brown and Ed Bueler and Constantine Khroulev and David Maxwell
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 "pism/util/Grid.hh"
20#include "pism/stressbalance/ssa/SSAFEM.hh"
21#include "pism/util/fem/Quadrature.hh"
22#include "pism/util/fem/DirichletData.hh"
23#include "pism/util/Mask.hh"
24#include "pism/basalstrength/basal_resistance.hh"
25#include "pism/rheology/FlowLaw.hh"
26#include "pism/util/pism_options.hh"
27#include "pism/util/error_handling.hh"
28#include "pism/util/Vars.hh"
29#include "pism/stressbalance/StressBalance.hh"
30#include "pism/geometry/Geometry.hh"
31
32#include "pism/util/node_types.hh"
33#include "pism/util/pism_utilities.hh" // average_water_column_pressure()
34#include "pism/util/Interpolation1D.hh"
35#include "pism/util/petscwrappers/DM.hh"
36#include "pism/util/petscwrappers/Vec.hh"
37#include "pism/util/petscwrappers/Viewer.hh"
38
39namespace pism {
40namespace stressbalance {
41
42/** The Q1 finite element SSA solver.
43 *
44 *
45 *
46 */
47SSAFEM::SSAFEM(std::shared_ptr<const Grid> grid)
48 : SSA(grid),
49 m_bc_mask(grid, "bc_mask"),
50 m_bc_values(grid, "_bc"),
51 m_gc(*m_config),
52 m_coefficients(grid, "ssa_coefficients", array::WITH_GHOSTS, 1),
53 m_node_type(m_grid, "node_type"),
54 m_boundary_integral(m_grid, "boundary_integral"),
55 m_element_index(*grid),
56 m_q1_element(*grid, fem::Q1Quadrature4()) {
59
60 const double ice_density = m_config->get_number("constants.ice.density");
61 m_alpha = 1 - ice_density / m_config->get_number("constants.sea_water.density");
62 m_rho_g = ice_density * m_config->get_number("constants.standard_gravity");
63
64 m_driving_stress_x = NULL;
65 m_driving_stress_y = NULL;
66
67 PetscErrorCode ierr;
68
69 m_dirichletScale = 1.0;
70 m_beta_ice_free_bedrock = m_config->get_number("basal_resistance.beta_ice_free_bedrock");
71
72 ierr = SNESCreate(m_grid->com, m_snes.rawptr());
73 PISM_CHK(ierr, "SNESCreate");
74
75 // Set the SNES callbacks to call into our compute_local_function and compute_local_jacobian.
77 m_callback_data.ssa = this;
78
79 auto dm = m_velocity_global.dm();
80
81 ierr = DMDASNESSetFunctionLocal(*dm, INSERT_VALUES,
82#if PETSC_VERSION_LT(3,21,0)
83 (DMDASNESFunction)function_callback,
84#else
85 (DMDASNESFunctionFn*)function_callback,
86#endif
88 PISM_CHK(ierr, "DMDASNESSetFunctionLocal");
89
90 ierr = DMDASNESSetJacobianLocal(*dm,
91#if PETSC_VERSION_LT(3,21,0)
92 (DMDASNESJacobian)jacobian_callback,
93#else
94 (DMDASNESJacobianFn*)jacobian_callback,
95#endif
97 PISM_CHK(ierr, "DMDASNESSetJacobianLocal");
98
99 ierr = DMSetMatType(*dm, "baij");
100 PISM_CHK(ierr, "DMSetMatType");
101
102 ierr = DMSetApplicationContext(*dm, &m_callback_data);
103 PISM_CHK(ierr, "DMSetApplicationContext");
104
105 ierr = SNESSetDM(m_snes, *dm);
106 PISM_CHK(ierr, "SNESSetDM");
107
108 // Default of maximum 200 iterations; possibly overridden by command line options
109 int snes_max_it = 200;
110 ierr = SNESSetTolerances(m_snes, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, snes_max_it,
111 PETSC_DEFAULT);
112 PISM_CHK(ierr, "SNESSetTolerances");
113
114 ierr = SNESSetFromOptions(m_snes);
115 PISM_CHK(ierr, "SNESSetFromOptions");
116
118 "node types: interior, boundary, exterior"); // no units or standard name
119
120 // Element::nodal_values() expects a ghosted array::Scalar. Ghosts if this field are
121 // never assigned to and not communicated, though.
123 "residual contribution from lateral boundaries"); // no units or standard name
124}
125
126// Initialize the solver, called once by the client before use.
128
130
131 // Use explicit driving stress if provided.
132 if (m_grid->variables().is_available("ssa_driving_stress_x") and
133 m_grid->variables().is_available("ssa_driving_stress_y")) {
134 m_driving_stress_x = m_grid->variables().get_2d_scalar("ssa_driving_stress_x");
135 m_driving_stress_y = m_grid->variables().get_2d_scalar("ssa_driving_stress_y");
136 }
137
138 m_log->message(2, " [using the SNES-based finite element method implementation]\n");
139
140 // process command-line options
141 {
142 m_dirichletScale = 1.0e9;
143 m_dirichletScale = options::Real(m_sys, "-ssa_fe_dirichlet_scale",
144 "Enforce Dirichlet conditions with this additional scaling",
145 "1", m_dirichletScale);
146 }
147
148 // On restart, SSA::init() reads the SSA velocity from a PISM output file
149 // into array::Vector "velocity". We use that field as an initial guess.
150 // If we are not restarting from a PISM file, "velocity" is identically zero,
151 // and the call below clears m_velocity_global.
152
154}
155
156/** Solve the SSA system of equations.
157
158 The FEM solver computes values of the various coefficients (most notably: the vertically-averaged
159 ice hardness) at each of the element nodes *only once*.
160
161 When running in an ice-model context, at each time step, SSA::update() is called, which calls
162 SSAFEM::solve(). Since coefficients have generally changed between time steps, we need to recompute
163 coefficients. On the other hand, in the context of inversion, coefficients will not change between
164 iteration and there is no need to recompute their values.
165
166 So there are two different solve methods, SSAFEM::solve() and SSAFEM::solve_nocache(). The only
167 difference is that SSAFEM::solve() recomputes the cached values of the coefficients before calling
168 SSAFEM::solve_nocache().
169 */
170void SSAFEM::solve(const Inputs &inputs) {
171
172 auto reason = solve_with_reason(inputs);
173 if (reason->failed()) {
175 "SSAFEM solve failed to converge (SNES reason %s)",
176 reason->description().c_str());
177 }
178
179 if (m_log->get_threshold() > 2) {
180 m_stdout_ssa += "SSAFEM converged (SNES reason " + reason->description() + ")";
181 }
182}
183
184std::shared_ptr<TerminationReason> SSAFEM::solve_with_reason(const Inputs &inputs) {
185
186 // Set up the system to solve.
187 cache_inputs(inputs);
188
189 return solve_nocache();
190}
191
192//! Solve the SSA without first recomputing the values of coefficients at quad
193//! points. See the disccusion of SSAFEM::solve for more discussion.
194std::shared_ptr<TerminationReason> SSAFEM::solve_nocache() {
195 PetscErrorCode ierr;
196
197 m_epsilon_ssa = m_config->get_number("stress_balance.ssa.epsilon");
198
199 options::String filename("-ssa_view", "");
200 if (filename.is_set()) {
201 petsc::Viewer viewer;
202 ierr = PetscViewerASCIIOpen(m_grid->com, filename->c_str(), viewer.rawptr());
203 PISM_CHK(ierr, "PetscViewerASCIIOpen");
204
205 ierr = PetscViewerASCIIPrintf(viewer, "SNES before SSASolve_FE\n");
206 PISM_CHK(ierr, "PetscViewerASCIIPrintf");
207
208 ierr = SNESView(m_snes, viewer);
209 PISM_CHK(ierr, "SNESView");
210
211 ierr = PetscViewerASCIIPrintf(viewer, "solution vector before SSASolve_FE\n");
212 PISM_CHK(ierr, "PetscViewerASCIIPrintf");
213
214 ierr = VecView(m_velocity_global.vec(), viewer);
215 PISM_CHK(ierr, "VecView");
216 }
217
218 m_stdout_ssa.clear();
219 if (m_log->get_threshold() >= 2) {
220 m_stdout_ssa = " SSA: ";
221 }
222
223 // Solve:
224 ierr = SNESSolve(m_snes, NULL, m_velocity_global.vec());
225 PISM_CHK(ierr, "SNESSolve");
226
227 // See if it worked.
228 SNESConvergedReason snes_reason;
229 ierr = SNESGetConvergedReason(m_snes, &snes_reason); PISM_CHK(ierr, "SNESGetConvergedReason");
230
231 std::shared_ptr<TerminationReason> reason(new SNESTerminationReason(snes_reason));
232 if (not reason->failed()) {
233
234 // Extract the solution back from SSAX to velocity and communicate.
236
237 bool view_solution = options::Bool("-ssa_view_solution", "view solution of the SSA system");
238 if (view_solution) {
239 petsc::Viewer viewer;
240 ierr = PetscViewerASCIIOpen(m_grid->com, filename->c_str(), viewer.rawptr());
241 PISM_CHK(ierr, "PetscViewerASCIIOpen");
242
243 ierr = PetscViewerASCIIPrintf(viewer, "solution vector after SSASolve\n");
244 PISM_CHK(ierr, "PetscViewerASCIIPrintf");
245
246 ierr = VecView(m_velocity_global.vec(), viewer);
247 PISM_CHK(ierr, "VecView");
248 }
249
250 }
251
252 return reason;
253}
254
255//! Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve.
256/**
257 This method is should be called after SSAFEM::init and whenever any geometry or temperature
258 related coefficients have changed. The method stores the values of the coefficients the nodes of
259 each element so that these are available to the residual and Jacobian evaluation methods.
260
261 We store the vertical average of the ice hardness to avoid re-doing this computation every
262 iteration.
263
264 In addition to coefficients at element nodes we store "node types" used to identify interior
265 elements, exterior elements, and boundary faces.
266*/
267void SSAFEM::cache_inputs(const Inputs &inputs) {
268
269 // Make copies of BC mask and BC values: they are needed in SNES callbacks and
270 // inputs.bc_{mask,values} are not available there.
271 if ((inputs.bc_mask != nullptr) and (inputs.bc_values != nullptr)) {
272 m_bc_mask.copy_from(*inputs.bc_mask);
274 } else {
275 m_bc_mask.set(0.0);
276 }
277
278 const std::vector<double> &z = m_grid->z();
279
281 inputs.enthalpy,
282 &inputs.geometry->ice_thickness,
283 &inputs.geometry->bed_elevation,
285 inputs.basal_yield_stress};
286
287 bool use_explicit_driving_stress = (m_driving_stress_x != NULL) && (m_driving_stress_y != NULL);
288 if (use_explicit_driving_stress) {
290 }
291
292 ParallelSection loop(m_grid->com);
293 try {
294 for (auto p = m_grid->points(); p; p.next()) {
295 const int i = p.i(), j = p.j();
296
297 double thickness = inputs.geometry->ice_thickness(i, j);
298
299 Vector2d tau_d;
300 if (use_explicit_driving_stress) {
301 tau_d.u = (*m_driving_stress_x)(i, j);
302 tau_d.v = (*m_driving_stress_y)(i, j);
303 } else {
304 // tau_d above is set to zero by the Vector2d
305 // constructor, but is not used
306 }
307
308 const double *enthalpy = inputs.enthalpy->get_column(i, j);
309 double hardness = rheology::averaged_hardness(*m_flow_law, thickness,
310 m_grid->kBelowHeight(thickness),
311 z.data(), enthalpy);
312
313 Coefficients c;
314 c.thickness = thickness;
315 c.bed = inputs.geometry->bed_elevation(i, j);
316 c.sea_level = inputs.geometry->sea_level_elevation(i, j);
317 c.tauc = (*inputs.basal_yield_stress)(i, j);
318 c.hardness = hardness;
319 c.driving_stress = tau_d;
320
321 m_coefficients(i, j) = c;
322 } // loop over owned grid points
323 } catch (...) {
324 loop.failed();
325 }
326 loop.check();
327
328 m_coefficients.update_ghosts();
329
330 const bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
331 if (use_cfbc) {
332 // Note: the call below uses ghosts of inputs.geometry->ice_thickness.
334 m_config->get_number("stress_balance.ice_free_thickness_standard"),
336 } else {
338 }
339
340 cache_residual_cfbc(inputs);
341
342}
343
344//! Compute quadrature point values of various coefficients given a quadrature `Q` and nodal values.
346 const Coefficients *x,
347 int *mask,
348 double *thickness,
349 double *tauc,
350 double *hardness) const {
351 // quadrature size
352 unsigned int n = E.n_pts();
353
354 for (unsigned int q = 0; q < n; q++) {
355 double
356 bed = 0.0,
357 sea_level = 0.0;
358
359 thickness[q] = 0.0;
360 tauc[q] = 0.0;
361 hardness[q] = 0.0;
362
363 for (int k = 0; k < E.n_chi(); k++) {
364 const fem::Germ &psi = E.chi(q, k);
365
366 thickness[q] += psi.val * x[k].thickness;
367 bed += psi.val * x[k].bed;
368 sea_level += psi.val * x[k].sea_level;
369 tauc[q] += psi.val * x[k].tauc;
370 hardness[q] += psi.val * x[k].hardness;
371 }
372
373 mask[q] = m_gc.mask(sea_level, bed, thickness[q]);
374 }
375}
376
377//! Compute gravitational driving stress at quadrature points.
378//! Uses explicitly-provided nodal values.
380 const Coefficients *x,
381 Vector2d *result) {
382 const unsigned int n = E.n_pts();
383
384 for (unsigned int q = 0; q < n; q++) {
385 result[q] = 0.0;
386
387 for (int k = 0; k < E.n_chi(); k++) {
388 const fem::Germ &psi = E.chi(q, k);
389 result[q] += psi.val * x[k].driving_stress;
390 }
391 }
392}
393
394/** Compute the the driving stress at quadrature points.
395 The surface elevation @f$ h @f$ is defined by
396 @f{equation}{
397 h =
398 \begin{cases}
399 b + H, & \mathrm{grounded},\\
400 z_{sl} + \alpha H, & \mathrm{shelf},
401 \end{cases}
402 @f}
403 where @f$ b @f$ is the bed elevation, @f$ H @f$ is the ice thickness, and @f$ z_{sl} @f$ is the
404 sea level, and @f$ \alpha @f$ is the ratio of densities of ice and sea water.
405
406 Note that if @f$ b @f$, @f$ H @f$, and @f$ z_{sl} @f$ are bilinear (or @f$ C^{1} @f$ in
407 general), @f$ h @f$ is continuous but not continuously differentiable. In
408 other words, the gradient of @f$ h @f$ is not continuous, so near the
409 grounding line we cannot compute the gravitational driving stress
410 @f$ \tau_{d} = - \rho g H \nabla h @f$ using the @f$ Q_1 @f$ or @f$ P_1 @f$ FE basis
411 expansion of @f$ h @f$.
412
413 We differentiate the equation above instead, getting
414 @f{equation}{
415 \nabla h =
416 \begin{cases}
417 \nabla b + \nabla H, & \mathrm{grounded},\\
418 \nabla z_{sl} + \alpha \nabla H, & \mathrm{shelf}.
419 \end{cases}
420 @f}
421
422 and use basis expansions of @f$ \nabla b @f$ and @f$ \nabla H @f$.
423
424 Overall, we have
425
426 @f{align*}{
427 \tau_{d} &= \rho g H \nabla h\\
428 &=
429 - \rho g H
430 \begin{cases}
431 \nabla b + \nabla H, & \mathrm{grounded},\\
432 \alpha \nabla H, & \mathrm{shelf},
433 \end{cases}
434 @f}
435
436 because @f$ z = z_{sl} @f$ defines the geoid surface and so *its gradient
437 does not contribute to the driving stress*.
438*/
440 const Coefficients *x,
441 Vector2d *result) const {
442 const unsigned int n = E.n_pts();
443
444 for (unsigned int q = 0; q < n; q++) {
445 double
446 sea_level = 0.0,
447 b = 0.0,
448 b_x = 0.0,
449 b_y = 0.0,
450 H = 0.0,
451 H_x = 0.0,
452 H_y = 0.0;
453
454 result[q] = 0.0;
455
456 for (int k = 0; k < E.n_chi(); k++) {
457 const fem::Germ &psi = E.chi(q, k);
458
459 b += psi.val * x[k].bed;
460 b_x += psi.dx * x[k].bed;
461 b_y += psi.dy * x[k].bed;
462
463 H += psi.val * x[k].thickness;
464 H_x += psi.dx * x[k].thickness;
465 H_y += psi.dy * x[k].thickness;
466
467 sea_level += psi.val * x[k].sea_level;
468 }
469
470 const int M = m_gc.mask(sea_level, b, H);
471 const bool grounded = mask::grounded(M);
472
473 const double
474 pressure = m_rho_g * H,
475 h_x = grounded ? b_x + H_x : m_alpha * H_x,
476 h_y = grounded ? b_y + H_y : m_alpha * H_y;
477
478 result[q].u = - pressure * h_x;
479 result[q].v = - pressure * h_y;
480 }
481}
482
483
484/** @brief Compute the "(regularized effective viscosity) x (ice thickness)" and effective viscous
485 * bed strength from the current solution, at a single quadrature point.
486 *
487 * @param[in] thickness ice thickness
488 * @param[in] hardness ice hardness
489 * @param[in] mask cell type mask
490 * @param[in] tauc basal yield stress
491 * @param[in] U the value of the solution
492 * @param[in] U_x x-derivatives of velocity components
493 * @param[in] U_y y-derivatives of velocity components
494 * @param[out] nuH product of the ice viscosity and thickness @f$ \nu H @f$
495 * @param[out] dnuH derivative of @f$ \nu H @f$ with respect to the
496 * second invariant @f$ \gamma @f$. Set to NULL if
497 * not desired.
498 * @param[out] beta basal drag coefficient @f$ \beta @f$
499 * @param[out] dbeta derivative of @f$ \beta @f$ with respect to the
500 * second invariant @f$ \gamma @f$. Set to NULL if
501 * not desired.
502 */
503void SSAFEM::PointwiseNuHAndBeta(double thickness,
504 double hardness,
505 int mask,
506 double tauc,
507 const Vector2d &U,
508 const Vector2d &U_x,
509 const Vector2d &U_y,
510 double *nuH, double *dnuH,
511 double *beta, double *dbeta) {
512
513 if (thickness < strength_extension->get_min_thickness()) {
515 if (dnuH != nullptr) {
516 *dnuH = 0;
517 }
518 } else {
519 m_flow_law->effective_viscosity(hardness, secondInvariant_2D(U_x, U_y),
520 nuH, dnuH);
521
522 *nuH = m_epsilon_ssa + *nuH * thickness;
523
524 if (dnuH != nullptr) {
525 *dnuH *= thickness;
526 }
527 }
528
529 if (mask::grounded_ice(mask)) {
530 m_basal_sliding_law->drag_with_derivative(tauc, U.u, U.v, beta, dbeta);
531 } else {
532 *beta = 0;
533
534 if (mask::ice_free_land(mask)) {
536 }
537
538 if (dbeta != nullptr) {
539 *dbeta = 0;
540 }
541 }
542}
543
544
545//! Compute and cache residual contributions from the integral over the lateral boundary.
546/**
547
548 This method computes FIXME.
549
550 */
552
553 // Compute 1D (linear) shape function values at quadrature points on element boundaries.
554 // These values do not depend on the element type or dimensions of a particular physical
555 // element.
556 //
557 // Note: the number of quadrature points (2) is hard-wired below.
558 const unsigned int Nq = 2;
559 double chi_b[Nq][Nq];
560 {
561 // interval length does not matter here
562 fem::Gaussian2 Q(1.0);
563
564 for (int i = 0; i < 2; ++i) { // 2 functions
565 for (int j = 0; j < 2; ++j) { // 2 quadrature points
566 chi_b[i][j] = fem::linear::chi(i, Q.point(j)).val;
567 }
568 }
569 }
570
571 const unsigned int Nk = fem::q1::n_chi;
572
573 using fem::P1Element2;
575 P1Element2 p1_element[Nk] = {P1Element2(*m_grid, Q_p1, 0),
576 P1Element2(*m_grid, Q_p1, 1),
577 P1Element2(*m_grid, Q_p1, 2),
578 P1Element2(*m_grid, Q_p1, 3)};
579
580 using mask::ocean;
581
582 const bool
583 use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
584
585 const double
586 ice_density = m_config->get_number("constants.ice.density"),
587 ocean_density = m_config->get_number("constants.sea_water.density"),
588 standard_gravity = m_config->get_number("constants.standard_gravity");
589
590 // Reset the boundary integral so that all values are overwritten.
592
593 if (not use_cfbc) {
594 // If CFBC is not used then we're done.
595 return;
596 }
597
599 &inputs.geometry->ice_thickness,
600 &inputs.geometry->bed_elevation,
603
604 // Iterate over the elements.
605 const int
606 xs = m_element_index.xs,
607 xm = m_element_index.xm,
608 ys = m_element_index.ys,
609 ym = m_element_index.ym;
610
611 ParallelSection loop(m_grid->com);
612 try {
613 for (int j = ys; j < ys + ym; j++) {
614 for (int i = xs; i < xs + xm; i++) {
615 // Initialize the map from global to element degrees of freedom.
616 m_q1_element.reset(i, j);
617
618 int node_type[Nk];
620
621 fem::Element2 *E = nullptr;
622 {
623 auto type = fem::element_type(node_type);
624
625 switch (type) {
627 continue;
628 case fem::ELEMENT_Q:
629 E = &m_q1_element;
630 break;
631 default:
632 E = &p1_element[type];
633
634 E->reset(i, j);
635
636 // get node types *again*, this time at the nodes of this P1 element
637 E->nodal_values(m_node_type, node_type);
638 }
639 }
640
641 // residual contributions at element nodes
642 std::vector<Vector2d> I(Nk);
643
644 double H_nodal[Nk];
645 E->nodal_values(inputs.geometry->ice_thickness.array(), H_nodal);
646
647 double b_nodal[Nk];
648 E->nodal_values(inputs.geometry->bed_elevation.array(), b_nodal);
649
650 double sl_nodal[Nk];
651 E->nodal_values(inputs.geometry->sea_level_elevation.array(), sl_nodal);
652
653 // storage for values of test functions on a side of the element
654 double psi[2] = {0.0, 0.0};
655
656 const unsigned int n_sides = E->n_sides();
657 // loop over element sides
658 for (unsigned int s = 0; s < n_sides; ++s) {
659
660 // initialize the quadrature and compute quadrature weights
662
663 // nodes incident to the current side
664 const int
665 n0 = s,
666 n1 = (n0 + 1) % n_sides;
667
668 if (not (node_type[n0] == NODE_BOUNDARY and
669 node_type[n1] == NODE_BOUNDARY)) {
670 // not a boundary side; skip it
671 continue;
672 }
673
674 for (unsigned int q = 0; q < Nq; ++q) {
675
676 double W = Q.weight(q);
677
678 // test functions at nodes incident to the current side, evaluated at the
679 // quadrature point q
680 psi[0] = chi_b[0][q];
681 psi[1] = chi_b[1][q];
682
683 // Compute ice thickness and bed elevation at a quadrature point. This uses a 1D basis
684 // expansion on the current side.
685 const double
686 H = H_nodal[n0] * psi[0] + H_nodal[n1] * psi[1],
687 bed = b_nodal[n0] * psi[0] + b_nodal[n1] * psi[1],
688 sea_level = sl_nodal[n0] * psi[0] + sl_nodal[n1] * psi[1];
689
690 // vertically averaged pressures at a quadrature point
691 double
692 P_ice = 0.5 * ice_density * standard_gravity * H,
693 P_water = average_water_column_pressure(H, bed, sea_level,
694 ice_density, ocean_density,
695 standard_gravity);
696
697 // vertical integral of the pressure difference
698 double dP = H * (P_ice - P_water);
699 // FIXME: implement melange pressure forcing
700
701 // This integral contributes to the residual at 2 nodes (the ones incident to the
702 // current side). This is is written in a way that allows *adding* (... += ...) the
703 // boundary contribution in the residual computation.
704 //
705 // FIXME: I need to include the special case corresponding to ice margins next
706 // to fjord walls, nunataks, etc. In this case dP == 0.
707 //
708 // FIXME: set pressure difference to zero at grounded locations at domain
709 // boundaries.
710 I[n0] += W * (- psi[0] * dP) * E->normal(s);
711 I[n1] += W * (- psi[1] * dP) * E->normal(s);
712 } // q-loop
713
714 } // loop over element sides
715
717
718 } // i-loop
719 } // j-loop
720 } catch (...) {
721 loop.failed();
722 }
723 loop.check();
724}
725
726
727//! Implements the callback for computing the residual.
728/*!
729 * Compute the residual \f[r_{ij}= G(x, \psi_{ij}) \f] where \f$G\f$
730 * is the weak form of the SSA, \f$x\f$ is the current approximate
731 * solution, and the \f$\psi_{ij}\f$ are test functions.
732 *
733 * The weak form of the SSA system is
734 */
735void SSAFEM::compute_local_function(Vector2d const *const *const velocity_global,
736 Vector2d **residual_global) {
737
738 const bool use_explicit_driving_stress = (m_driving_stress_x != NULL) && (m_driving_stress_y != NULL);
739
740 const bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
741
742 const unsigned int Nk = fem::q1::n_chi;
743 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
744
745 using fem::P1Element2;
747 P1Element2 p1_element[Nk] = {P1Element2(*m_grid, Q_p1, 0),
748 P1Element2(*m_grid, Q_p1, 1),
749 P1Element2(*m_grid, Q_p1, 2),
750 P1Element2(*m_grid, Q_p1, 3)};
751
753
754 // Set the boundary contribution of the residual. This is computed at the nodes, so we don't want
755 // to set it using Element::add_contribution() because that would lead to
756 // double-counting. Also note that without CFBC m_boundary_integral is exactly zero.
757 for (auto p = m_grid->points(); p; p.next()) {
758 const int i = p.i(), j = p.j();
759
760 residual_global[j][i] = m_boundary_integral(i, j);
761 }
762
763 // Start access to Dirichlet data if present.
765
766 // Storage for the current solution and its derivatives at quadrature points.
767 Vector2d U[Nq_max], U_x[Nq_max], U_y[Nq_max];
768
769 // Iterate over the elements.
770 const int
771 xs = m_element_index.xs,
772 xm = m_element_index.xm,
773 ys = m_element_index.ys,
774 ym = m_element_index.ym;
775
776 ParallelSection loop(m_grid->com);
777 try {
778 for (int j = ys; j < ys + ym; j++) {
779 for (int i = xs; i < xs + xm; i++) {
780
781 fem::Element2 *E = nullptr;
782 {
783 m_q1_element.reset(i, j);
784
785 int node_type[Nk];
787
788 auto type = fem::element_type(node_type);
789
790 if (use_cfbc) {
791 if (type == fem::ELEMENT_EXTERIOR) {
792 // skip exterior elements
793 continue;
794 }
795
796 if (type == fem::ELEMENT_Q) {
797 E = &m_q1_element;
798 } else {
799 E = &p1_element[type];
800
801 E->reset(i, j);
802 }
803 } else {
804 // if use_cfbc == false all elements are interior and Q1
805 E = &m_q1_element;
806 }
807 }
808
809 // Number of quadrature points.
810 const unsigned int Nq = E->n_pts();
811
812 // Storage for the solution and residuals at element nodes.
813 Vector2d residual[Nk];
814
815 int mask[Nq_max];
816 double thickness[Nq_max];
817 double tauc[Nq_max];
818 double hardness[Nq_max];
819 Vector2d tau_d[Nq_max];
820
821 {
822 Coefficients coeffs[Nk];
823 E->nodal_values(m_coefficients.array(), coeffs);
824
825 quad_point_values(*E, coeffs, mask, thickness, tauc, hardness);
826
827 if (use_explicit_driving_stress) {
828 explicit_driving_stress(*E, coeffs, tau_d);
829 } else {
830 driving_stress(*E, coeffs, tau_d);
831 }
832 }
833
834 {
835 // Obtain the value of the solution at the nodes
836 Vector2d velocity_nodal[Nk];
837 E->nodal_values(velocity_global, velocity_nodal);
838
839 // These values now need to be adjusted if some nodes in the element have Dirichlet data.
840 if (dirichlet_data) {
841 // Set elements of velocity_nodal that correspond to Dirichlet nodes to prescribed
842 // values.
843 dirichlet_data.enforce(*E, velocity_nodal);
844 // mark Dirichlet nodes in E so that they are not touched by
845 // add_contribution() below
846 dirichlet_data.constrain(*E);
847 }
848
849 // Compute the solution values and its gradient at the quadrature points.
850 E->evaluate(velocity_nodal, // input
851 U, U_x, U_y); // outputs
852 }
853
854 // Zero out the element-local residual in preparation for updating it.
855 for (unsigned int k = 0; k < Nk; k++) {
856 residual[k] = 0;
857 }
858
859 // loop over quadrature points:
860 for (unsigned int q = 0; q < Nq; q++) {
861
862 auto W = E->weight(q);
863
864 double eta = 0.0, beta = 0.0;
865 PointwiseNuHAndBeta(thickness[q], hardness[q], mask[q], tauc[q],
866 U[q], U_x[q], U_y[q], // inputs
867 &eta, NULL, &beta, NULL); // outputs
868
869 // The next few lines compute the actual residual for the element.
870 const Vector2d tau_b = U[q] * (- beta); // basal shear stress
871
872 const double
873 u_x = U_x[q].u,
874 v_y = U_y[q].v,
875 u_y_plus_v_x = U_y[q].u + U_x[q].v;
876
877 // Loop over test functions.
878 for (int k = 0; k < E->n_chi(); k++) {
879 const fem::Germ &psi = E->chi(q, k);
880
881 residual[k].u += W * (eta * (psi.dx * (4.0 * u_x + 2.0 * v_y) + psi.dy * u_y_plus_v_x)
882 - psi.val * (tau_b.u + tau_d[q].u));
883 residual[k].v += W * (eta * (psi.dx * u_y_plus_v_x + psi.dy * (2.0 * u_x + 4.0 * v_y))
884 - psi.val * (tau_b.v + tau_d[q].v));
885 } // k (test functions)
886 } // q (quadrature points)
887
888 E->add_contribution(residual, residual_global);
889 } // i-loop
890 } // j-loop
891 } catch (...) {
892 loop.failed();
893 }
894 loop.check();
895
896 // Until now we have not touched rows in the residual corresponding to Dirichlet data.
897 // We fix this now.
898 if (dirichlet_data) {
899 dirichlet_data.fix_residual(velocity_global, residual_global);
900 }
901
902 if (use_cfbc) {
903 // Prescribe homogeneous Dirichlet B.C. at ice-free nodes. This uses the fact that m_node_type
904 // can be used as a "Dirichlet B.C. mask", i.e. ice-free nodes (and only ice-free nodes) are
905 // marked with ones.
906 fem::DirichletData_Vector dirichlet_ice_free(&m_node_type, NULL, m_dirichletScale);
907 dirichlet_ice_free.fix_residual_homogeneous(residual_global);
908 }
909
910 monitor_function(velocity_global, residual_global);
911}
912
913void SSAFEM::monitor_function(Vector2d const *const *const velocity_global,
914 Vector2d const *const *const residual_global) {
915 PetscErrorCode ierr;
916 bool monitorFunction = options::Bool("-ssa_monitor_function", "monitor the SSA residual");
917 if (not monitorFunction) {
918 return;
919 }
920
921 ierr = PetscPrintf(m_grid->com,
922 "SSA Solution and Function values (pointwise residuals)\n");
923 PISM_CHK(ierr, "PetscPrintf");
924
925 ParallelSection loop(m_grid->com);
926 try {
927 for (auto p = m_grid->points(); p; p.next()) {
928 const int i = p.i(), j = p.j();
929
930 ierr = PetscSynchronizedPrintf(m_grid->com,
931 "[%2d, %2d] u=(%12.10e, %12.10e) f=(%12.4e, %12.4e)\n",
932 i, j,
933 velocity_global[j][i].u, velocity_global[j][i].v,
934 residual_global[j][i].u, residual_global[j][i].v);
935 PISM_CHK(ierr, "PetscSynchronizedPrintf");
936 }
937 } catch (...) {
938 loop.failed();
939 }
940 loop.check();
941
942 ierr = PetscSynchronizedFlush(m_grid->com, NULL);
943 PISM_CHK(ierr, "PetscSynchronizedFlush");
944}
945
946
947//! Implements the callback for computing the Jacobian.
948/*!
949 Compute the Jacobian
950
951 @f[ J_{ij}{kl} \frac{d r_{ij}}{d x_{kl}}= G(x, \psi_{ij}) @f]
952
953 where \f$G\f$ is the weak form of the SSA, \f$x\f$ is the current
954 approximate solution, and the \f$\psi_{ij}\f$ are test functions.
955
956*/
957void SSAFEM::compute_local_jacobian(Vector2d const *const *const velocity_global, Mat Jac) {
958
959 const bool use_cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
960
961 const unsigned int Nk = fem::q1::n_chi;
962 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
963
964 using fem::P1Element2;
966 P1Element2 p1_element[Nk] = {P1Element2(*m_grid, Q_p1, 0),
967 P1Element2(*m_grid, Q_p1, 1),
968 P1Element2(*m_grid, Q_p1, 2),
969 P1Element2(*m_grid, Q_p1, 3)};
970
971 // Zero out the Jacobian in preparation for updating it.
972 PetscErrorCode ierr = MatZeroEntries(Jac);
973 PISM_CHK(ierr, "MatZeroEntries");
974
976
977 // Start access to Dirichlet data if present.
979
980 // Storage for the current solution at quadrature points.
981 Vector2d U[Nq_max], U_x[Nq_max], U_y[Nq_max];
982
983 // Loop through all the elements.
984 int
985 xs = m_element_index.xs,
986 xm = m_element_index.xm,
987 ys = m_element_index.ys,
988 ym = m_element_index.ym;
989
990 ParallelSection loop(m_grid->com);
991 try {
992 for (int j = ys; j < ys + ym; j++) {
993 for (int i = xs; i < xs + xm; i++) {
994
995 fem::Element2 *E = nullptr;
996 {
997 m_q1_element.reset(i, j);
998
999 int node_type[Nk];
1001
1002 auto type = fem::element_type(node_type);
1003
1004 if (use_cfbc) {
1005 if (type == fem::ELEMENT_EXTERIOR) {
1006 // skip exterior elements
1007 continue;
1008 }
1009
1010 if (type == fem::ELEMENT_Q) {
1011 E = &m_q1_element;
1012 } else {
1013 E = &p1_element[type];
1014
1015 E->reset(i, j);
1016 }
1017 } else {
1018 // if use_cfbc == false all elements are interior and Q1
1019 E = &m_q1_element;
1020 }
1021 }
1022
1023 // Number of quadrature points.
1024 const unsigned int
1025 Nq = E->n_pts(),
1026 n_chi = E->n_chi();
1027
1028 int mask[Nq_max];
1029 double thickness[Nq_max];
1030 double tauc[Nq_max];
1031 double hardness[Nq_max];
1032
1033 {
1034 Coefficients coeffs[Nk];
1035 E->nodal_values(m_coefficients.array(), coeffs);
1036
1037 quad_point_values(*E, coeffs,
1038 mask, thickness, tauc, hardness);
1039 }
1040
1041 {
1042 // Values of the solution at the nodes of the current element.
1043 Vector2d velocity_nodal[Nk];
1044 E->nodal_values(velocity_global, velocity_nodal);
1045
1046 // These values now need to be adjusted if some nodes in the element have
1047 // Dirichlet data.
1048 if (dirichlet_data) {
1049 dirichlet_data.enforce(*E, velocity_nodal);
1050 dirichlet_data.constrain(*E);
1051 }
1052 // Compute the values of the solution at the quadrature points.
1053 E->evaluate(velocity_nodal, U, U_x, U_y);
1054 }
1055
1056 // Element-local Jacobian matrix (there are Nk vector valued degrees
1057 // of freedom per element, for a total of (2*Nk)*(2*Nk) = 64
1058 // entries in the local Jacobian.
1059 double K[2*Nk][2*Nk];
1060 // Build the element-local Jacobian.
1061 ierr = PetscMemzero(K, sizeof(K));
1062 PISM_CHK(ierr, "PetscMemzero");
1063
1064 for (unsigned int q = 0; q < Nq; q++) {
1065
1066 const double
1067 W = E->weight(q),
1068 u = U[q].u,
1069 v = U[q].v,
1070 u_x = U_x[q].u,
1071 v_y = U_y[q].v,
1072 u_y_plus_v_x = U_y[q].u + U_x[q].v;
1073
1074 double eta = 0.0, deta = 0.0, beta = 0.0, dbeta = 0.0;
1075 PointwiseNuHAndBeta(thickness[q], hardness[q], mask[q], tauc[q],
1076 U[q], U_x[q], U_y[q],
1077 &eta, &deta, &beta, &dbeta);
1078
1079 for (unsigned int l = 0; l < n_chi; l++) { // Trial functions
1080
1081 // Current trial function and its derivatives:
1082 const fem::Germ &phi = E->chi(q, l);
1083
1084 // Derivatives of \gamma with respect to u_l and v_l:
1085 const double
1086 gamma_u = (2.0 * u_x + v_y) * phi.dx + 0.5 * u_y_plus_v_x * phi.dy,
1087 gamma_v = 0.5 * u_y_plus_v_x * phi.dx + (u_x + 2.0 * v_y) * phi.dy;
1088
1089 // Derivatives of \eta = \nu*H with respect to u_l and v_l:
1090 const double
1091 eta_u = deta * gamma_u,
1092 eta_v = deta * gamma_v;
1093
1094 // Derivatives of the basal shear stress term (\tau_b):
1095 const double
1096 taub_xu = -dbeta * u * u * phi.val - beta * phi.val, // x-component, derivative with respect to u_l
1097 taub_xv = -dbeta * u * v * phi.val, // x-component, derivative with respect to u_l
1098 taub_yu = -dbeta * v * u * phi.val, // y-component, derivative with respect to v_l
1099 taub_yv = -dbeta * v * v * phi.val - beta * phi.val; // y-component, derivative with respect to v_l
1100
1101 for (unsigned int k = 0; k < n_chi; k++) { // Test functions
1102
1103 // Current test function and its derivatives:
1104
1105 const fem::Germ &psi = E->chi(q, k);
1106
1107 if (eta == 0) {
1108 ierr = PetscPrintf(PETSC_COMM_SELF, "eta=0 i %d j %d q %d k %d\n", i, j, q, k);
1109 PISM_CHK(ierr, "PetscPrintf");
1110 }
1111
1112 // u-u coupling
1113 K[k*2 + 0][l*2 + 0] += W * (eta_u * (psi.dx * (4 * u_x + 2 * v_y) + psi.dy * u_y_plus_v_x)
1114 + eta * (4 * psi.dx * phi.dx + psi.dy * phi.dy) - psi.val * taub_xu);
1115 // u-v coupling
1116 K[k*2 + 0][l*2 + 1] += W * (eta_v * (psi.dx * (4 * u_x + 2 * v_y) + psi.dy * u_y_plus_v_x)
1117 + eta * (2 * psi.dx * phi.dy + psi.dy * phi.dx) - psi.val * taub_xv);
1118 // v-u coupling
1119 K[k*2 + 1][l*2 + 0] += W * (eta_u * (psi.dx * u_y_plus_v_x + psi.dy * (2 * u_x + 4 * v_y))
1120 + eta * (psi.dx * phi.dy + 2 * psi.dy * phi.dx) - psi.val * taub_yu);
1121 // v-v coupling
1122 K[k*2 + 1][l*2 + 1] += W * (eta_v * (psi.dx * u_y_plus_v_x + psi.dy * (2 * u_x + 4 * v_y))
1123 + eta * (psi.dx * phi.dx + 4 * psi.dy * phi.dy) - psi.val * taub_yv);
1124
1125 } // l
1126 } // k
1127 } // q
1128 E->add_contribution(&K[0][0], Jac);
1129 } // j
1130 } // i
1131 } catch (...) {
1132 loop.failed();
1133 }
1134 loop.check();
1135
1136
1137 // Until now, the rows and columns correspoinding to Dirichlet data
1138 // have not been set. We now put an identity block in for these
1139 // unknowns. Note that because we have takes steps to not touching
1140 // these columns previously, the symmetry of the Jacobian matrix is
1141 // preserved.
1142 if (dirichlet_data) {
1143 dirichlet_data.fix_jacobian(Jac);
1144 }
1145
1146 if (use_cfbc) {
1147 // Prescribe homogeneous Dirichlet B.C. at ice-free nodes. This uses the fact that m_node_type
1148 // can be used as a "Dirichlet B.C. mask", i.e. ice-free nodes (and only ice-free nodes) are
1149 // marked with ones.
1150 fem::DirichletData_Vector dirichlet_ice_free(&m_node_type, NULL, m_dirichletScale);
1151 dirichlet_ice_free.fix_jacobian(Jac);
1152 }
1153
1154 ierr = MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
1155 PISM_CHK(ierr, "MatAssemblyBegin");
1156
1157 ierr = MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
1158 PISM_CHK(ierr, "MatAssemblyEnd");
1159
1160 ierr = MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1161 PISM_CHK(ierr, "MatSetOption");
1162
1163 ierr = MatSetOption(Jac, MAT_SYMMETRIC, PETSC_TRUE);
1164 PISM_CHK(ierr, "MatSetOption");
1165
1166 monitor_jacobian(Jac);
1167}
1168
1170 PetscErrorCode ierr;
1171 bool mon_jac = options::Bool("-ssa_monitor_jacobian", "monitor the SSA Jacobian");
1172
1173 if (not mon_jac) {
1174 return;
1175 }
1176
1177 // iter has to be a PetscInt because it is used in the
1178 // SNESGetIterationNumber() call below.
1179 PetscInt iter = 0;
1180 ierr = SNESGetIterationNumber(m_snes, &iter);
1181 PISM_CHK(ierr, "SNESGetIterationNumber");
1182
1183 auto file_name = pism::printf("PISM_SSAFEM_J%d.m", (int)iter);
1184
1185 m_log->message(2,
1186 "writing Matlab-readable file for SSAFEM system A xsoln = rhs to file `%s' ...\n",
1187 file_name.c_str());
1188
1189 petsc::Viewer viewer(m_grid->com);
1190
1191 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);
1192 PISM_CHK(ierr, "PetscViewerSetType");
1193
1194 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1195 PISM_CHK(ierr, "PetscViewerPushFormat");
1196
1197 ierr = PetscViewerFileSetName(viewer, file_name.c_str());
1198 PISM_CHK(ierr, "PetscViewerFileSetName");
1199
1200 ierr = PetscObjectSetName((PetscObject) Jac, "A");
1201 PISM_CHK(ierr, "PetscObjectSetName");
1202
1203 ierr = MatView(Jac, viewer);
1204 PISM_CHK(ierr, "MatView");
1205
1206 ierr = PetscViewerPopFormat(viewer);
1207 PISM_CHK(ierr, "PetscViewerPopFormat");
1208}
1209
1210//!
1211PetscErrorCode SSAFEM::function_callback(DMDALocalInfo *info,
1212 Vector2d const *const *const velocity,
1213 Vector2d **residual,
1214 CallbackData *fe) {
1215 try {
1216 (void) info;
1217 fe->ssa->compute_local_function(velocity, residual);
1218 } catch (...) {
1219 MPI_Comm com = MPI_COMM_SELF;
1220 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)fe->da, &com); CHKERRQ(ierr);
1222 SETERRQ(com, 1, "A PISM callback failed");
1223 }
1224 return 0;
1225}
1226
1227PetscErrorCode SSAFEM::jacobian_callback(DMDALocalInfo *info,
1228 Vector2d const *const *const velocity,
1229 Mat A, Mat J, CallbackData *fe) {
1230 try {
1231 (void) A;
1232 (void) info;
1234 } catch (...) {
1235 MPI_Comm com = MPI_COMM_SELF;
1236 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)fe->da, &com); CHKERRQ(ierr);
1238 SETERRQ(com, 1, "A PISM callback failed");
1239 }
1240 return 0;
1241}
1242
1243} // end of namespace stressbalance
1244} // end of namespace pism
const units::System::Ptr m_sys
unit system used by this component
Definition Component.hh:160
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
int mask(double sea_level, double bed, double thickness) const
Definition Mask.hh:138
array::Scalar1 sea_level_elevation
Definition Geometry.hh:48
array::Scalar2 ice_thickness
Definition Geometry.hh:51
array::Scalar2 bed_elevation
Definition Geometry.hh:47
virtual void drag_with_derivative(double tauc, double vx, double vy, double *drag, double *ddrag) const
Compute the drag coefficient and its derivative with respect to .
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
VariableMetadata & long_name(const std::string &input)
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
void copy_from(const Array2D< T > &source)
Definition Array2D.hh:73
void add(double alpha, const Array2D< T > &x)
Definition Array2D.hh:65
double * get_column(int i, int j)
Definition Array3D.cc:121
void set_interpolation_type(InterpolationType type)
Definition Array.cc:178
petsc::Vec & vec() const
Definition Array.cc:310
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
std::shared_ptr< petsc::DM > dm() const
Definition Array.cc:324
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
void fix_residual(Vector2d const *const *x_global, Vector2d **r_global)
void enforce(const Element2 &element, Vector2d *x_e)
void fix_residual_homogeneous(Vector2d **r)
void constrain(Element2 &element)
Constrain element, i.e. ensure that quadratures do not contribute to Dirichlet nodes by marking corre...
void evaluate(const T *x, T *vals, T *dx, T *dy)
Given nodal values, compute the values and partial derivatives at the quadrature points.
Definition Element.hh:195
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
double side_length(unsigned int side) const
Definition Element.hh:185
void add_contribution(const T *local, T **y_global) const
Add the values of element-local contributions y to the global vector y_global.
Definition Element.hh:238
void nodal_values(const array::Scalar &x_global, int *result) const
Get nodal values of an integer mask.
Definition Element.cc:185
unsigned int n_sides() const
Definition Element.hh:181
Vector2d normal(unsigned int side) const
Definition Element.hh:176
int xm
total number of elements to loop over in the x-direction.
int ym
total number of elements to loop over in the y-direction.
int ys
y-coordinate of the first element to loop over.
int xs
x-coordinate of the first element to loop over.
int n_chi() const
Definition Element.hh:65
double weight(unsigned int q) const
Weight of the quadrature point q
Definition Element.hh:85
const Germ & chi(unsigned int q, unsigned int k) const
Definition Element.hh:73
unsigned int n_pts() const
Number of quadrature points.
Definition Element.hh:80
The mapping from global to local degrees of freedom.
Definition Element.hh:61
P1 element embedded in Q1Element2.
Definition Element.hh:273
double weight(int k) const
Definition Quadrature.hh:72
QuadPoint point(int k) const
Definition Quadrature.hh:68
bool is_set() const
Definition options.hh:35
const array::Scalar * basal_yield_stress
const array::Scalar * bc_mask
const array::Array3D * enthalpy
const array::Vector * bc_values
std::shared_ptr< TerminationReason > solve_with_reason(const Inputs &inputs)
Definition SSAFEM.cc:184
void PointwiseNuHAndBeta(double thickness, double hardness, int mask, double tauc, const Vector2d &U, const Vector2d &U_x, const Vector2d &U_y, double *nuH, double *dnuH, double *beta, double *dbeta)
Compute the "(regularized effective viscosity) x (ice thickness)" and effective viscous bed strength ...
Definition SSAFEM.cc:503
void driving_stress(const fem::Element &E, const Coefficients *x, Vector2d *result) const
Definition SSAFEM.cc:439
const array::Scalar * m_driving_stress_x
Definition SSAFEM.hh:142
void cache_inputs(const Inputs &inputs)
Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve.
Definition SSAFEM.cc:267
static PetscErrorCode function_callback(DMDALocalInfo *info, Vector2d const *const *velocity, Vector2d **residual, CallbackData *fe)
SNES callbacks.
Definition SSAFEM.cc:1211
SSAFEM(std::shared_ptr< const Grid > g)
Definition SSAFEM.cc:47
static void explicit_driving_stress(const fem::Element &E, const Coefficients *x, Vector2d *result)
Definition SSAFEM.cc:379
CallbackData m_callback_data
Definition SSAFEM.hh:122
void compute_local_jacobian(Vector2d const *const *velocity, Mat J)
Implements the callback for computing the Jacobian.
Definition SSAFEM.cc:957
array::Vector1 m_boundary_integral
Boundary integral (CFBC contribution to the residual).
Definition SSAFEM.hh:129
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, Vector2d const *const *velocity, Mat A, Mat J, CallbackData *fe)
Definition SSAFEM.cc:1227
void monitor_jacobian(Mat Jac)
Definition SSAFEM.cc:1169
void monitor_function(Vector2d const *const *velocity_global, Vector2d const *const *residual_global)
Definition SSAFEM.cc:913
array::Scalar1 m_bc_mask
Definition SSAFEM.hh:67
virtual void solve(const Inputs &inputs)
Definition SSAFEM.cc:170
fem::Q1Element2 m_q1_element
Definition SSAFEM.hh:136
GeometryCalculator m_gc
Definition SSAFEM.hh:70
virtual void init_impl()
Initialize a generic regular-grid SSA solver.
Definition SSAFEM.cc:127
std::shared_ptr< TerminationReason > solve_nocache()
Definition SSAFEM.cc:194
void cache_residual_cfbc(const Inputs &inputs)
Compute and cache residual contributions from the integral over the lateral boundary.
Definition SSAFEM.cc:551
void compute_local_function(Vector2d const *const *velocity, Vector2d **residual)
Implements the callback for computing the residual.
Definition SSAFEM.cc:735
fem::ElementIterator m_element_index
Definition SSAFEM.hh:135
void quad_point_values(const fem::Element &E, const Coefficients *x, int *mask, double *thickness, double *tauc, double *hardness) const
Compute quadrature point values of various coefficients given a quadrature Q and nodal values.
Definition SSAFEM.cc:345
const array::Scalar * m_driving_stress_y
Definition SSAFEM.hh:143
array::Scalar1 m_node_type
Storage for node types (interior, boundary, exterior).
Definition SSAFEM.hh:127
array::Vector1 m_bc_values
Definition SSAFEM.hh:68
array::Array2D< Coefficients > m_coefficients
Definition SSAFEM.hh:74
double get_notional_strength() const
Returns strength = (viscosity times thickness).
Definition SSA.cc:59
array::Vector m_velocity_global
Definition SSA.hh:133
SSAStrengthExtension * strength_extension
Definition SSA.hh:115
std::string m_stdout_ssa
Definition SSA.hh:131
virtual void init_impl()
Initialize a generic regular-grid SSA solver.
Definition SSA.cc:101
PISM's SSA solver.
Definition SSA.hh:110
const array::Vector1 & velocity() const
Get the thickness-advective 2D velocity.
std::shared_ptr< rheology::FlowLaw > m_flow_law
IceBasalResistancePlasticLaw * m_basal_sliding_law
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
const double phi
Definition exactTestK.c:41
#define n
Definition exactTestM.c:37
Germ chi(unsigned int k, const QuadPoint &pt)
Linear basis functions on the interval [-1, -1].
Definition FEM.cc:30
const int n_chi
Definition FEM.hh:191
const unsigned int MAX_QUADRATURE_SIZE
Definition FEM.hh:173
ElementType element_type(const int node_type[q1::n_chi])
Definition FEM.cc:106
@ ELEMENT_Q
Definition FEM.hh:203
@ ELEMENT_EXTERIOR
Definition FEM.hh:205
bool grounded_ice(int M)
Definition Mask.hh:51
bool ice_free_land(int M)
Definition Mask.hh:64
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition Mask.hh:44
bool ocean(int M)
An ocean cell (floating ice or ice-free).
Definition Mask.hh:40
bool Bool(const std::string &option, const std::string &description)
Definition options.cc:190
double averaged_hardness(const FlowLaw &ice, double thickness, unsigned int kbelowH, const double *zlevels, const double *enthalpy)
Computes vertical average of B(E, p) ice hardness, namely .
Definition FlowLaw.cc:213
@ NODE_BOUNDARY
Definition node_types.hh:35
@ NODE_INTERIOR
Definition node_types.hh:34
static double secondInvariant_2D(const Vector2d &U_x, const Vector2d &U_y)
Definition FlowLaw.hh:45
void compute_node_types(const array::Scalar1 &ice_thickness, double thickness_threshold, array::Scalar &result)
Definition node_types.cc:62
double average_water_column_pressure(double ice_thickness, double bed, double floatation_level, double rho_ice, double rho_water, double g)
std::string printf(const char *format,...)
static const double k
Definition exactTestP.cc:42
void handle_fatal_errors(MPI_Comm com)
double dy
Function derivative with respect to y.
Definition FEM.hh:157
double val
Function value.
Definition FEM.hh:153
double dx
Function deriviative with respect to x.
Definition FEM.hh:155
Struct for gathering the value and derivative of a function at a point.
Definition FEM.hh:151
Adaptor for gluing SNESDAFormFunction callbacks to an SSAFEM.
Definition SSAFEM.hh:116
double tauc
basal yield stress
Definition SSAFEM.hh:60
Vector2d driving_stress
prescribed gravitational driving stress
Definition SSAFEM.hh:64