PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
jacobian.cc
Go to the documentation of this file.
1 /* Copyright (C) 2020, 2021 PISM Authors
2  *
3  * This file is part of PISM.
4  *
5  * PISM is free software; you can redistribute it and/or modify it under the
6  * terms of the GNU General Public License as published by the Free Software
7  * Foundation; either version 3 of the License, or (at your option) any later
8  * version.
9  *
10  * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13  * details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with PISM; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #include "Blatter.hh"
21 
22 #include "pism/basalstrength/basal_resistance.hh"
23 #include "pism/rheology/FlowLaw.hh"
24 #include "pism/util/node_types.hh"
25 
26 #include "util/DataAccess.hh"
27 #include "util/grid_hierarchy.hh" // grid_transpose(), grid_z()
28 
29 namespace pism {
30 namespace stressbalance {
31 
32 /*!
33  * Computes the Jacobian contribution of the "main" part of the Blatter system.
34  */
36  const Vector2 *u_nodal,
37  const double *B_nodal,
38  double K[16][16]) {
39  int Nk = fem::q13d::n_chi;
40 
41  Vector2
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  element.evaluate(u_nodal, u, u_x, u_y, u_z);
50  element.evaluate(B_nodal, B);
51 
52  // loop over all quadrature points
53  for (int q = 0; q < element.n_pts(); ++q) {
54  auto W = element.weight(q) / m_scaling;
55 
56  double
57  ux = u_x[q].u,
58  uy = u_y[q].u,
59  uz = u_z[q].u,
60  vx = u_x[q].v,
61  vy = u_y[q].v,
62  vz = u_z[q].v;
63 
64  double gamma = (ux * ux + vy * vy + ux * vy +
65  0.25 * ((uy + vx) * (uy + vx) + uz * uz + vz * vz));
66 
67  double eta, deta;
68  m_flow_law->effective_viscosity(B[q], gamma, m_viscosity_eps, &eta, &deta);
69 
70  // add the enhancement factor
71  eta *= m_E_viscosity;
72  deta *= m_E_viscosity;
73 
74  // loop over test and trial functions, computing the upper-triangular part of
75  // the element Jacobian
76  for (int t = 0; t < Nk; ++t) {
77  auto psi = element.chi(q, t);
78  for (int s = t; s < Nk; ++s) {
79  auto phi = element.chi(q, s);
80 
81  // partial derivatives of gamma with respect to u_i and v_i
82  double
83  gamma_u = 2.0 * ux * phi.dx + vy * phi.dx + 0.5 * phi.dy * (uy + vx) + 0.5 * uz * phi.dz,
84  gamma_v = 2.0 * vy * phi.dy + ux * phi.dy + 0.5 * phi.dx * (uy + vx) + 0.5 * vz * phi.dz;
85 
86  // partial derivatives of eta with respect to u_i and v_i, using chain rule
87  double
88  eta_u = deta * gamma_u,
89  eta_v = deta * gamma_v;
90 
91  // F_u = grad(psi) . (4ux + 2vy, uy + vx, uz) and
92  // F_v = grad(psi) . (uy + vx, 4vy + 2ux, vz)
93  double
94  F_u = (psi.dx * (4.0 * ux + 2.0 * vy) + psi.dy * (uy + vx) + psi.dz * uz),
95  F_v = (psi.dx * (uy + vx) + psi.dy * (4.0 * vy + 2.0 * ux) + psi.dz * vz);
96 
97  // partial derivatives of F_u with respect to u_i and v_i
98  double
99  F_uu = 4.0 * psi.dx * phi.dx + psi.dy * phi.dy + psi.dz * phi.dz,
100  F_uv = 2.0 * psi.dx * phi.dy + psi.dy * phi.dx;
101 
102  // partial derivatives of F_v with respect to u_i and v_i
103  double
104  F_vu = 2.0 * psi.dy * phi.dx + psi.dx * phi.dy,
105  F_vv = 4.0 * psi.dy * phi.dy + psi.dx * phi.dx + psi.dz * phi.dz;
106 
107  K[t * 2 + 0][s * 2 + 0] += W * (eta * F_uu + eta_u * F_u);
108  K[t * 2 + 0][s * 2 + 1] += W * (eta * F_uv + eta_v * F_u);
109  K[t * 2 + 1][s * 2 + 0] += W * (eta * F_vu + eta_u * F_v);
110  K[t * 2 + 1][s * 2 + 1] += W * (eta * F_vv + eta_v * F_v);
111  }
112  }
113  } // end of the loop over q
114 
115 }
116 
117 /*!
118  * Compute the Jacobian contribution of the basal boundary condition.
119  *
120  * This method implements basal sliding.
121  */
123  const double *tauc_nodal,
124  const double *f_nodal,
125  const Vector2 *u_nodal,
126  double K[16][16]) {
127  int Nk = fem::q13d::n_chi;
128 
129  Vector2 *u = m_work2[0];
130 
131  double
132  *tauc = m_work[0],
133  *floatation = m_work[1];
134 
135  face.evaluate(u_nodal, u);
136  face.evaluate(tauc_nodal, tauc);
137  face.evaluate(f_nodal, floatation);
138 
139  for (int q = 0; q < face.n_pts(); ++q) {
140  auto W = face.weight(q) / m_scaling;
141 
142  bool grounded = floatation[q] <= 0.0;
143  double beta = 0.0, dbeta = 0.0;
144  if (grounded) {
145  m_basal_sliding_law->drag_with_derivative(tauc[q], u[q].u, u[q].v, &beta, &dbeta);
146  }
147 
148  // loop over all test functions
149  for (int t = 0; t < Nk; ++t) {
150  auto psi = face.chi(q, t);
151  for (int s = 0; s < Nk; ++s) {
152  auto phi = face.chi(q, s);
153 
154  double p = psi * phi;
155 
156  K[t * 2 + 0][s * 2 + 0] += W * p * (beta + dbeta * u[q].u * u[q].u);
157  K[t * 2 + 0][s * 2 + 1] += W * p * dbeta * u[q].u * u[q].v;
158  K[t * 2 + 1][s * 2 + 0] += W * p * dbeta * u[q].v * u[q].u;
159  K[t * 2 + 1][s * 2 + 1] += W * p * (beta + dbeta * u[q].v * u[q].v);
160  }
161  }
162  }
163 }
164 
165 /*!
166  * Set the Jacobian to identity at Dirichlet nodes.
167  */
168 void Blatter::jacobian_dirichlet(const DMDALocalInfo &info, Parameters **P, Mat J) {
169  PetscErrorCode ierr;
170 
171  // Dirichlet scaling
172  Vector2 scaling = {1.0, 1.0};
173 
174  // take care of Dirichlet nodes (both explicit and grid points outside the domain)
175  //
176  // here we loop over all the *owned* nodes
177  for (int j = info.ys; j < info.ys + info.ym; j++) {
178  for (int i = info.xs; i < info.xs + info.xm; i++) {
179  for (int k = info.zs; k < info.zs + info.zm; k++) {
180  if ((int)P[j][i].node_type == NODE_EXTERIOR or dirichlet_node(info, {i, j, k})) {
181 
182  double identity[4] = {scaling.u, 0, 0, scaling.v};
183 
184  MatStencil row;
185  row.i = k; // STORAGE_ORDER
186  row.j = i; // STORAGE_ORDER
187  row.k = j; // STORAGE_ORDER
188  ierr = MatSetValuesBlockedStencil(J, 1, &row, 1, &row, identity, ADD_VALUES);
189  PISM_CHK(ierr, "MatSetValuesBlockedStencil"); // this may throw
190  }
191  }
192  }
193  }
194 }
195 
196 /*!
197  * Compute the Jacobian matrix.
198  */
199 void Blatter::compute_jacobian(DMDALocalInfo *petsc_info,
200  const Vector2 ***X, Mat A, Mat J) {
201  auto info = grid_transpose(*petsc_info);
202 
203  // Zero out the Jacobian in preparation for updating it.
204  PetscErrorCode ierr = MatZeroEntries(J);
205  PISM_CHK(ierr, "MatZeroEntries");
206 
207  ierr = MatSetOption(A, MAT_SUBSET_OFF_PROC_ENTRIES, PETSC_TRUE);
208  PISM_CHK(ierr, "MatSetOption");
209 
210  ierr = MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
211  PISM_CHK(ierr, "MatSetOption");
212 
213  ierr = MatSetOption(J, MAT_SYMMETRIC, PETSC_TRUE);
214  PISM_CHK(ierr, "MatSetOption");
215 
216  ierr = PetscObjectSetName((PetscObject)J, "bp_jacobian");
217  PISM_CHK(ierr, "PetscObjectSetName");
218 
219  // Stencil width of 1 is not very important, but if info.sw > 1 will lead to more
220  // redundant computation (we would be looping over elements that don't contribute to any
221  // owned nodes).
222  assert(info.sw == 1);
223 
224  // horizontal grid spacing is the same on all multigrid levels
225  double
226  x_min = m_grid->x0() - m_grid->Lx(),
227  y_min = m_grid->y0() - m_grid->Ly(),
228  dx = m_grid->dx(),
229  dy = m_grid->dy();
230 
231  fem::Q1Element3 element(info,
233  dx, dy, x_min, y_min);
234 
235  // Maximum number of nodes per element
236  const int Nk = fem::q13d::n_chi;
237  assert(element.n_chi() <= Nk);
238  assert(element.n_pts() <= m_Nq);
239 
240  // scalar quantities
241  double z[Nk];
242  double floatation[Nk], bottom_elevation[Nk], ice_thickness[Nk];
243  double B_nodal[Nk], basal_yield_stress[Nk];
244  int node_type[Nk];
245 
246  // 2D vector quantities
247  Vector2 velocity[Nk];
248 
249  // note: we use info.da below because ice hardness is on the grid corresponding to the
250  // current multigrid level
251  //
252  // FIXME: This communicates ghosts of ice hardness
253  DataAccess<double***> hardness(info.da, 3, GHOSTED);
254 
256  auto *P = m_parameters.array();
257 
258  // loop over all the elements that have at least one owned node
259  for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
260  for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
261 
262  // Initialize 2D geometric info at element nodes
263  nodal_parameter_values(element, P, i, j,
264  node_type,
265  bottom_elevation,
266  ice_thickness,
267  NULL,
268  NULL);
269 
270  // skip ice-free (exterior) columns
271  if (exterior_element(node_type)) {
272  continue;
273  }
274 
275  for (int k = info.gzs; k < info.gzs + info.gzm - 1; k++) {
276 
277  // Element-local Jacobian matrix (there are Nk vector valued degrees of freedom
278  // per element, for a total of Nk*Nk = 64 entries in the local Jacobian.
279  double K[2*Nk][2*Nk];
280  memset(K, 0, sizeof(K));
281 
282  // Compute coordinates of the nodes of this element and fetch node types.
283  for (int n = 0; n < Nk; ++n) {
284  auto I = element.local_to_global(i, j, k, n);
285 
286  z[n] = grid_z(bottom_elevation[n], ice_thickness[n], info.mz, I.k);
287  }
288 
289  // compute values of chi, chi_x, chi_y, chi_z and quadrature weights at quadrature
290  // points on this physical element
291  element.reset(i, j, k, z);
292 
293  // Get nodal values of ice velocity.
294  {
295  element.nodal_values(X, velocity);
296 
297  // Don't contribute to Dirichlet nodes
298  for (int n = 0; n < Nk; ++n) {
299  auto I = element.local_to_global(n);
300  if (dirichlet_node(info, I)) {
301  element.mark_row_invalid(n);
302  element.mark_col_invalid(n);
303  velocity[n] = u_bc(element.x(n), element.y(n), element.z(n));
304  }
305  }
306  }
307 
308  element.nodal_values((double***)hardness, B_nodal);
309 
310  jacobian_f(element, velocity, B_nodal, K);
311 
312  // basal boundary
313  if (k == 0) {
314  for (int n = 0; n < Nk; ++n) {
315  auto I = element.local_to_global(n);
316 
317  basal_yield_stress[n] = P[I.j][I.i].tauc;
318  floatation[n] = P[I.j][I.i].floatation;
319  }
320 
321  fem::Q1Element3Face *face = grounding_line(floatation) ? &m_face100 : &m_face4;
322 
323  face->reset(fem::q13d::FACE_BOTTOM, z);
324 
325  jacobian_basal(*face, basal_yield_stress, floatation, velocity, K);
326  }
327 
328  // fill the lower-triangular part of the element Jacobian using the fact that J is
329  // symmetric
330  for (int t = 0; t < Nk; ++t) {
331  for (int s = 0; s < t; ++s) {
332  K[t * 2 + 0][s * 2 + 0] = K[s * 2 + 0][t * 2 + 0];
333  K[t * 2 + 1][s * 2 + 0] = K[s * 2 + 0][t * 2 + 1];
334  K[t * 2 + 0][s * 2 + 1] = K[s * 2 + 1][t * 2 + 0];
335  K[t * 2 + 1][s * 2 + 1] = K[s * 2 + 1][t * 2 + 1];
336  }
337  }
338 
339  element.add_contribution(&K[0][0], J);
340  } // end of the loop over i
341  } // end of the loop over j
342  } // end of the loop over k
343 
344  jacobian_dirichlet(info, P, J);
345 
346  ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); PISM_CHK(ierr, "MatAssemblyBegin");
347  ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); PISM_CHK(ierr, "MatAssemblyEnd");
348  if (A != J) {
349  ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); PISM_CHK(ierr, "MatAssemblyBegin");
350  ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); PISM_CHK(ierr, "MatAssemblyEnd");
351  }
352 }
353 
354 PetscErrorCode Blatter::jacobian_callback(DMDALocalInfo *info,
355  const Vector2 ***x,
356  Mat A, Mat J,
357  Blatter *solver) {
358  try {
359  solver->compute_jacobian(info, x, A, J);
360  } catch (...) {
361  MPI_Comm com = solver->grid()->com;
362  handle_fatal_errors(com);
363  SETERRQ(com, 1, "A PISM callback failed");
364  }
365  return 0;
366 }
367 
368 } // end of namespace stressbalance
369 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
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 .
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 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
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
void mark_col_invalid(unsigned int k)
Mark that the column corresponding to local degree of freedom k should not be updated when inserting ...
Definition: Element.cc:218
3D Q1 element
Definition: Element.hh:351
void jacobian_dirichlet(const DMDALocalInfo &info, Parameters **P, Mat J)
Definition: jacobian.cc:168
virtual void jacobian_basal(const fem::Q1Element3Face &face, const double *tauc_nodal, const double *f_nodal, const Vector2 *u_nodal, double K[2 *fem::q13d::n_chi][2 *fem::q13d::n_chi])
Definition: jacobian.cc:122
void compute_jacobian(DMDALocalInfo *info, const Vector2 ***x, Mat A, Mat J)
Definition: jacobian.cc:199
Vector2 m_work2[m_n_work][m_Nq]
Definition: Blatter.hh:110
fem::Q1Element3Face m_face4
Definition: Blatter.hh:112
virtual void jacobian_f(const fem::Q1Element3 &element, const Vector2 *u_nodal, const double *B_nodal, double K[2 *fem::q13d::n_chi][2 *fem::q13d::n_chi])
Definition: jacobian.cc:35
double m_work[m_n_work][m_Nq]
Definition: Blatter.hh:108
static const int m_Nq
Definition: Blatter.hh:105
fem::Q1Element3Face m_face100
Definition: Blatter.hh:113
virtual Vector2 u_bc(double x, double y, double z) const
Definition: Blatter.cc:131
static bool exterior_element(const int *node_type)
Definition: Blatter.cc:145
static bool grounding_line(const double *F)
Definition: Blatter.cc:165
virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex &I)
Definition: Blatter.cc:123
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, const Vector2 ***x, Mat A, Mat J, Blatter *solver)
Definition: jacobian.cc:354
IceModelVec2< Parameters > m_parameters
Definition: Blatter.hh:79
virtual void nodal_parameter_values(const fem::Q1Element3 &element, Parameters **P, int i, int j, int *node_type, double *bottom, double *thickness, double *surface, double *sea_level) const
Definition: Blatter.cc:627
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)
const double phi
Definition: exactTestK.c:41
#define n
Definition: exactTestM.c:37
const int n_chi
Number of shape functions on a Q1 element.
Definition: FEM.hh:213
static Vector3 row(const double A[3][3], size_t k)
Definition: Element.cc:46
static double K(double psi_x, double psi_y, double speed, double epsilon)
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:43
double grid_z(double b, double H, int Mz, int k)
@ 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:45
void handle_fatal_errors(MPI_Comm com)
const int J[]
Definition: ssafd_code.cc:34
const int I[]
Definition: ssafd_code.cc:24