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