PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
residual.cc
Go to the documentation of this file.
1 /* Copyright (C) 2020, 2021, 2023 PISM Authors
2  *
3  * This file is part of PISM.
4  *
5  * PISM is free software; you can redistribute it and/or modify it under the
6  * terms of the GNU General Public License as published by the Free Software
7  * Foundation; either version 3 of the License, or (at your option) any later
8  * version.
9  *
10  * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13  * details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with PISM; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #include "pism/stressbalance/blatter/Blatter.hh"
21 
22 #include "pism/basalstrength/basal_resistance.hh"
23 #include "pism/rheology/FlowLaw.hh"
24 #include "pism/util/node_types.hh"
25 #include "pism/stressbalance/blatter/util/DataAccess.hh"
26 #include "pism/stressbalance/blatter/util/grid_hierarchy.hh" // grid_transpose(), grid_z()
27 
28 namespace pism {
29 namespace stressbalance {
30 
31 static const Vector2d u_exterior = {0.0, 0.0};
32 
33 /*!
34  * Computes the residual contribution of the "main" part of the Blatter system.
35  */
37  const Vector2d *u_nodal,
38  const double *B_nodal,
39  Vector2d *residual) {
40 
41  Vector2d
42  *u = m_work2[0],
43  *u_x = m_work2[1],
44  *u_y = m_work2[2],
45  *u_z = m_work2[3];
46 
47  double *B = m_work[0];
48 
49  // evaluate u and its partial derivatives at quadrature points
50  element.evaluate(u_nodal, u, u_x, u_y, u_z);
51 
52  // evaluate B (ice hardness) at quadrature points
53  element.evaluate(B_nodal, B);
54 
55  // loop over all quadrature points
56  for (int q = 0; q < element.n_pts(); ++q) {
57  auto W = element.weight(q) / m_scaling;
58 
59  double
60  ux = u_x[q].u,
61  uy = u_y[q].u,
62  uz = u_z[q].u,
63  vx = u_x[q].v,
64  vy = u_y[q].v,
65  vz = u_z[q].v;
66 
67  double gamma = (ux * ux + vy * vy + ux * vy +
68  0.25 * ((uy + vx) * (uy + vx) + uz * uz + vz * vz));
69 
70  double eta;
71  m_flow_law->effective_viscosity(B[q], gamma, m_viscosity_eps, &eta, nullptr);
72 
73  // add the enhancement factor
74  eta *= m_E_viscosity;
75 
76  // loop over all test functions
77  for (int t = 0; t < element.n_chi(); ++t) {
78  const auto &psi = element.chi(q, t);
79 
80  residual[t].u += W * (eta * (psi.dx * (4.0 * ux + 2.0 * vy) +
81  psi.dy * (uy + vx) +
82  psi.dz * uz));
83  residual[t].v += W * (eta * (psi.dx * (uy + vx) +
84  psi.dy * (2.0 * ux + 4.0 * vy) +
85  psi.dz * vz));
86  }
87  }
88 }
89 
90 /*! Computes the residual contribution of the "source term".
91  *
92  * This term contains the driving stress.
93  */
95  const double *surface,
96  const double *bed,
97  Vector2d *residual) {
98 
99  // compute s_x and s_y
100  double
101  *s_x = m_work[0],
102  *s_y = m_work[1];
103 
104  if (m_eta_transform) {
105  // use the same storage for eta_x and eta_y
106  double
107  *eta_nodal = m_work[2],
108  *eta = m_work[3],
109  *eta_x = s_x,
110  *eta_y = s_y,
111  *eta_z = m_work[4];
112  double
113  *b = m_work[5],
114  *b_x = m_work[6],
115  *b_y = m_work[7],
116  *b_z = m_work[8];
117 
118  double
119  n = m_glen_exponent,
120  p = (2.0 * n + 2.0) / n;
121 
122  for (int k = 0; k < element.n_chi(); ++k) {
123  eta_nodal[k] = pow(surface[k] - bed[k], p);
124  }
125 
126  element.evaluate(eta_nodal, eta, eta_x, eta_y, eta_z);
127  element.evaluate(bed, b, b_x, b_y, b_z);
128 
129  for (int q = 0; q < element.n_pts(); ++q) {
130  double C = pow(eta[q], 1.0 / p - 1.0) / p;
131 
132  s_x[q] = C * eta_x[q] + b_x[q];
133  s_y[q] = C * eta_y[q] + b_y[q];
134  }
135  } else {
136  // these arrays are needed by the call below (but results are discarded)
137  double
138  *s = m_work[2],
139  *s_z = m_work[3];
140 
141  element.evaluate(surface, s, s_x, s_y, s_z);
142  }
143 
144  for (int q = 0; q < element.n_pts(); ++q) {
145  auto W = element.weight(q) / m_scaling;
146 
147  auto F = m_rho_ice_g * Vector2d(s_x[q], s_y[q]);
148 
149  // loop over all test functions
150  for (int t = 0; t < element.n_chi(); ++t) {
151  const auto &psi = element.chi(q, t);
152 
153  residual[t] += W * psi.val * F;
154  }
155  }
156 }
157 
158 /*!
159  * Computes the residual contribution of the basal boundary condition.
160  *
161  * This takes care of basal sliding.
162  */
164  const fem::Q1Element3Face &face,
165  const double *tauc_nodal,
166  const double *f_nodal,
167  const Vector2d *u_nodal,
168  Vector2d *residual) {
169 
170  Vector2d *u = m_work2[0];
171 
172  double
173  *tauc = m_work[0],
174  *floatation = m_work[1];
175 
176  face.evaluate(u_nodal, u);
177  face.evaluate(tauc_nodal, tauc);
178  face.evaluate(f_nodal, floatation);
179 
180  for (int q = 0; q < face.n_pts(); ++q) {
181  auto W = face.weight(q) / m_scaling;
182 
183  bool grounded = floatation[q] <= 0.0;
184  double beta = grounded ? m_basal_sliding_law->drag(tauc[q], u[q].u, u[q].v) : 0.0;
185 
186  // loop over all test functions
187  for (int t = 0; t < element.n_chi(); ++t) {
188  auto psi = face.chi(q, t);
189 
190  residual[t] += W * psi * beta * u[q];
191  }
192  }
193 }
194 
195 /*!
196  * Top surface contribution to the residual.
197  *
198  * Used by verification tests ONLY.
199  */
201  const fem::Q1Element3Face &face,
202  Vector2d *residual) {
203  (void) element;
204  (void) face;
205  (void) residual;
206  // In normal circumstances the top surface contribution is zero (natural BCs apply).
207 }
208 
209 
210 /*!
211  * Computes the residual contribution of lateral boundary conditions.
212  *
213  * This takes care of "calving front" stress boundary conditions.
214  *
215  * FIXME: make p_ocean an input from a parameterization of the melange back pressure.
216  */
218  const fem::Q1Element3Face &face,
219  const double *surface_nodal,
220  const double *z_nodal,
221  const double *sl_nodal,
222  Vector2d *residual) {
223  double
224  *z = m_work[0],
225  *s = m_work[1],
226  *sea_level = m_work[2];
227 
228  // compute physical coordinates of quadrature points on this face
229  face.evaluate(surface_nodal, s);
230  face.evaluate(z_nodal, z);
231  face.evaluate(sl_nodal, sea_level);
232 
233  // loop over all quadrature points
234  for (int q = 0; q < face.n_pts(); ++q) {
235  auto W = face.weight(q) / m_scaling;
236  auto N3 = face.normal(q);
237  Vector2d N = {N3.x, N3.y};
238 
239  double
240  ice_depth = s[q] - z[q],
241  p_ice = m_rho_ice_g * ice_depth,
242  water_depth = std::max(sea_level[q] - z[q], 0.0),
243  p_ocean = m_rho_ocean_g * water_depth;
244 
245  // loop over all test functions
246  for (int t = 0; t < element.n_chi(); ++t) {
247  auto psi = face.chi(q, t);
248 
249  residual[t] += - W * psi * (p_ice - p_ocean) * N;
250  }
251  }
252 }
253 
254 /*! Set the residual at Dirichlet locations
255  *
256  * Compute the residual at Dirichlet locations and reset the residual to zero elsewhere.
257  *
258  * Setting it to zero is necessary because we call DMDASNESSetFunctionLocal() with
259  * INSERT_VALUES.
260  *
261  */
262 void Blatter::residual_dirichlet(const DMDALocalInfo &info,
263  Parameters **P,
264  const Vector2d ***x,
265  Vector2d ***R) {
266 
267  double
268  x_min = m_grid->x0() - m_grid->Lx(),
269  y_min = m_grid->y0() - m_grid->Ly(),
270  dx = m_grid->dx(),
271  dy = m_grid->dy(),
272  scaling = 1.0;
273 
274  // Compute the residual at Dirichlet BC nodes and reset the residual to zero elsewhere.
275  //
276  // here we loop over all the *owned* nodes
277  for (int j = info.ys; j < info.ys + info.ym; j++) {
278  for (int i = info.xs; i < info.xs + info.xm; i++) {
279 
280  // compute the residual at ice-free map plane locations
281  if ((int)P[j][i].node_type == NODE_EXTERIOR) {
282  for (int k = info.zs; k < info.zs + info.zm; k++) {
283  R[j][i][k] = scaling * (x[j][i][k] - u_exterior); // STORAGE_ORDER
284  }
285  continue;
286  }
287 
288  for (int k = info.zs; k < info.zs + info.zm; k++) {
289  // reset to zero
290  R[j][i][k] = 0.0; // STORAGE_ORDER
291 
292  // compute the residual at Dirichlet nodes in verification tests
293  if (dirichlet_node(info, {i, j, k})) {
294  double
295  xx = x_min + i * dx,
296  yy = y_min + j * dy,
297  b = P[j][i].bed,
298  H = P[j][i].thickness,
299  zz = grid_z(b, H, info.mz, k);
300  Vector2d U_bc = u_bc(xx, yy, zz);
301  R[j][i][k] = scaling * (x[j][i][k] - U_bc); // STORAGE_ORDER
302  }
303  }
304  }
305  }
306 }
307 
308 /*!
309  * Computes the residual.
310  */
311 void Blatter::compute_residual(DMDALocalInfo *petsc_info,
312  const Vector2d ***X, Vector2d ***R) {
313  auto info = grid_transpose(*petsc_info);
314 
315  // Stencil width of 1 is not very important, but if info.sw > 1 will lead to more
316  // redundant computation (we would be looping over elements that don't contribute to any
317  // owned nodes).
318  assert(info.sw == 1);
319 
320  // horizontal grid spacing is the same on all multigrid levels
321  double
322  x_min = m_grid->x0() - m_grid->Lx(),
323  y_min = m_grid->y0() - m_grid->Ly(),
324  dx = m_grid->dx(),
325  dy = m_grid->dy();
326 
327  fem::Q1Element3 element(info, fem::Q13DQuadrature8(),
328  dx, dy, x_min, y_min);
329 
330  // Number of nodes per element.
331  const int Nk = fem::q13d::n_chi;
332  assert(element.n_chi() <= Nk);
333  assert(element.n_pts() <= m_Nq);
334 
335  // scalar quantities
336  double z[Nk];
337  double floatation[Nk], sea_level[Nk], bottom_elevation[Nk], ice_thickness[Nk], surface_elevation[Nk];
338  double B[Nk], basal_yield_stress[Nk];
339  int node_type[Nk];
340 
341  // 2D vector quantities
342  Vector2d velocity[Nk], R_nodal[Nk];
343 
344  // note: we use info.da below because ice hardness is on the grid corresponding to the
345  // current multigrid level
346  //
347  // FIXME: This communicates ghosts of ice hardness
348  DataAccess<double***> ice_hardness(info.da, 3, GHOSTED);
349 
351  auto *P = m_parameters.array();
352 
353  // Compute the residual at Dirichlet nodes and set it to zero elsewhere.
354  residual_dirichlet(info, P, X, R);
355 
356  // loop over all the elements that have at least one owned node
357  for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
358  for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
359 
360  // Initialize 2D geometric info at element nodes
361  nodal_parameter_values(element, P, i, j,
362  node_type,
363  bottom_elevation,
364  ice_thickness,
365  surface_elevation,
366  sea_level);
367 
368  // skip ice-free (exterior) elements
369  if (exterior_element(node_type)) {
370  continue;
371  }
372 
373  // loop over elements in a column
374  for (int k = info.gzs; k < info.gzs + info.gzm - 1; k++) {
375 
376  // Reset element residual to zero in preparation.
377  for (int n = 0; n < Nk; ++n) {
378  R_nodal[n] = 0.0;
379  }
380 
381  // Compute z coordinates of the nodes of this element
382  for (int n = 0; n < Nk; ++n) {
383  auto I = element.local_to_global(i, j, k, n);
384  z[n] = grid_z(bottom_elevation[n], ice_thickness[n], info.mz, I.k);
385  }
386 
387  // Compute values of chi, chi_x, chi_y, chi_z and quadrature weights at quadrature
388  // points on this physical element
389  element.reset(i, j, k, z);
390 
391  // Get nodal values of ice velocity.
392  {
393  element.nodal_values(X, velocity);
394 
395  // Take care of Dirichlet BC: don't contribute to Dirichlet nodes and set nodal
396  // values of the current iterate to Dirichlet BC values.
397  for (int n = 0; n < Nk; ++n) {
398  auto I = element.local_to_global(n);
399  if (dirichlet_node(info, I)) {
400  element.mark_row_invalid(n);
401  velocity[n] = u_bc(element.x(n), element.y(n), element.z(n));
402  }
403  }
404  }
405 
406  // Get nodal values of ice hardness.
407  element.nodal_values((double***)ice_hardness, B);
408 
409  // "main" part of the residual
410  residual_f(element, velocity, B, R_nodal);
411 
412  // the "source term" (driving stress)
413  residual_source_term(element, surface_elevation, bottom_elevation, R_nodal);
414 
415  // basal boundary
416  if (k == 0) {
417  for (int n = 0; n < Nk; ++n) {
418  auto I = element.local_to_global(n);
419 
420  basal_yield_stress[n] = P[I.j][I.i].tauc;
421  floatation[n] = P[I.j][I.i].floatation;
422  }
423 
424  // use an N*N-point equally-spaced quadrature at grounding lines
425  fem::Q1Element3Face *face = grounding_line(floatation) ? &m_face100 : &m_face4;
426  face->reset(fem::q13d::FACE_BOTTOM, z);
427 
428  residual_basal(element, *face, basal_yield_stress, floatation, velocity, R_nodal);
429  }
430 
431  // lateral boundary
432  // loop over all vertical faces (see fem::q13d::incident_nodes for the order)
433  for (int f = 0; f < 4; ++f) {
434  if (marine_boundary(f, node_type, bottom_elevation, sea_level)) {
435  // use an N*N-point equally-spaced quadrature for partially-submerged faces
436  fem::Q1Element3Face *face = (partially_submerged_face(f, z, sea_level) ?
437  &m_face100 : &m_face4);
438  face->reset(f, z);
439 
440  residual_lateral(element, *face, surface_elevation, z, sea_level, R_nodal);
441  }
442  } // end of the loop over element faces
443 
444  // top boundary (verification tests only)
445  if (k == info.mz - 2) {
447 
448  residual_surface(element, m_face4, R_nodal);
449  }
450 
451  element.add_contribution(R_nodal, R);
452  } // end of the loop over i
453  } // end of the loop over j
454  } // end of the loop over k
455 }
456 
457 PetscErrorCode Blatter::function_callback(DMDALocalInfo *info,
458  const Vector2d ***x, Vector2d ***f,
459  Blatter *solver) {
460  try {
461  solver->compute_residual(info, x, f);
462  } catch (...) {
463  MPI_Comm com = solver->grid()->com;
464  handle_fatal_errors(com);
465  SETERRQ(com, 1, "A PISM callback failed");
466  }
467  return 0;
468 }
469 
470 } // end of namespace stressbalance
471 } // end of namespace pism
std::shared_ptr< const Grid > grid() const
Definition: Component.cc:105
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
virtual double drag(double tauc, double vx, double vy) const
Compute the drag coefficient for the basal shear stress.
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2d.hh:29
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void nodal_values(T const *const *const *x_global, T *result) const
Extract nodal values for the element (i,j,k) from global array x_global into the element-local array ...
Definition: Element.hh:314
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:328
GlobalIndex local_to_global(int i, int j, int k, unsigned int n) const
Definition: Element.hh:298
void evaluate(const T *x, T *vals, T *dx, T *dy, T *dz) const
Given nodal values, compute the values and partial derivatives at the quadrature points.
Definition: Element.hh:278
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
double weight(unsigned int q) const
Definition: Element.hh:407
const double & chi(unsigned int q, unsigned int k) const
Definition: Element.hh:412
void evaluate(const T *x, T *result) const
Definition: Element.hh:427
const Vector3 & normal(unsigned int q) const
Definition: Element.hh:419
int n_pts() const
Number of quadrature points.
Definition: Element.hh:404
void reset(int face, const double *z)
Definition: Element.cc:534
void reset(int i, int j, int k, const double *z)
Definition: Element.cc:451
double y(int n) const
Definition: Element.hh:369
double z(int n) const
Definition: Element.hh:374
void mark_row_invalid(unsigned int k)
Mark that the row corresponding to local degree of freedom k should not be updated when inserting int...
Definition: Element.cc:212
double x(int n) const
Definition: Element.hh:364
3D Q1 element
Definition: Element.hh:351
virtual void residual_f(const fem::Q1Element3 &element, const Vector2d *u_nodal, const double *B_nodal, Vector2d *residual)
Definition: residual.cc:36
virtual Vector2d u_bc(double x, double y, double z) const
Definition: Blatter.cc:132
virtual void residual_source_term(const fem::Q1Element3 &element, const double *surface, const double *bed, Vector2d *residual)
Definition: residual.cc:94
virtual void residual_basal(const fem::Q1Element3 &element, const fem::Q1Element3Face &face, const double *tauc_nodal, const double *f_nodal, const Vector2d *u_nodal, Vector2d *residual)
Definition: residual.cc:163
void compute_residual(DMDALocalInfo *info, const Vector2d ***X, Vector2d ***R)
Definition: residual.cc:311
virtual bool marine_boundary(int face, const int *node_type, const double *ice_bottom, const double *sea_level)
Definition: Blatter.cc:225
fem::Q1Element3Face m_face4
Definition: Blatter.hh:112
static bool partially_submerged_face(int face, const double *z, const double *sea_level)
Definition: Blatter.cc:192
virtual void residual_lateral(const fem::Q1Element3 &element, const fem::Q1Element3Face &face, const double *surface_nodal, const double *z_nodal, const double *sl_nodal, Vector2d *residual)
Definition: residual.cc:217
double m_work[m_n_work][m_Nq]
Definition: Blatter.hh:108
void residual_dirichlet(const DMDALocalInfo &info, Parameters **P, const Vector2d ***x, Vector2d ***R)
Definition: residual.cc:262
static PetscErrorCode function_callback(DMDALocalInfo *info, const Vector2d ***x, Vector2d ***f, Blatter *solver)
Definition: residual.cc:457
static const int m_Nq
Definition: Blatter.hh:105
fem::Q1Element3Face m_face100
Definition: Blatter.hh:113
static bool exterior_element(const int *node_type)
Definition: Blatter.cc:146
static bool grounding_line(const double *F)
Definition: Blatter.cc:166
array::Array2D< Parameters > m_parameters
Definition: Blatter.hh:79
virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex &I)
Definition: Blatter.cc:124
virtual void residual_surface(const fem::Q1Element3 &element, const fem::Q1Element3Face &face, Vector2d *residual)
Definition: residual.cc:200
Vector2d m_work2[m_n_work][m_Nq]
Definition: Blatter.hh:110
virtual void nodal_parameter_values(const fem::Q1Element3 &element, Parameters **P, int i, int j, int *node_type, double *bottom, double *thickness, double *surface, double *sea_level) const
Definition: Blatter.cc:627
const array::Vector1 & velocity() const
Get the thickness-advective 2D velocity.
std::shared_ptr< rheology::FlowLaw > m_flow_law
IceBasalResistancePlasticLaw * m_basal_sliding_law
#define n
Definition: exactTestM.c:37
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
const int n_chi
Number of shape functions on a Q1 element.
Definition: FEM.hh:213
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:44
static const Vector2d u_exterior
Definition: residual.cc:31
double grid_z(double b, double H, int Mz, int k)
static double F(double SL, double B, double H, double alpha)
@ NODE_EXTERIOR
Definition: node_types.hh:36
@ GHOSTED
Definition: DataAccess.hh:30
DMDALocalInfo grid_transpose(const DMDALocalInfo &input)
static const double k
Definition: exactTestP.cc:42
void handle_fatal_errors(MPI_Comm com)
const int I[]
Definition: ssafd_code.cc:24
const int C[]
Definition: ssafd_code.cc:44