PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
IP_SSAHardavForwardProblem.cc
Go to the documentation of this file.
1 // Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 2020, 2021 David Maxwell and Constantine Khroulev
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 
20 #include "pism/basalstrength/basal_resistance.hh"
21 #include "pism/util/IceGrid.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/Vars.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/pism_utilities.hh"
27 #include "pism/rheology/FlowLaw.hh"
28 #include "pism/geometry/Geometry.hh"
29 #include "pism/stressbalance/StressBalance.hh"
30 #include "pism/util/petscwrappers/DM.hh"
31 #include "pism/util/petscwrappers/Vec.hh"
32 
33 namespace pism {
34 namespace inverse {
35 
38  : SSAFEM(g),
39  m_stencil_width(1),
40  m_zeta(NULL),
41  m_dzeta_local(m_grid, "d_zeta_local", WITH_GHOSTS, m_stencil_width),
42  m_fixed_design_locations(NULL),
43  m_design_param(tp),
44  m_du_global(m_grid, "linearization work vector (sans ghosts)",
45  WITHOUT_GHOSTS, m_stencil_width),
46  m_du_local(m_grid, "linearization work vector (with ghosts)",
47  WITH_GHOSTS, m_stencil_width),
48  m_hardav(m_grid, "hardav", WITH_GHOSTS, m_stencil_width),
49  m_element_index(*m_grid),
50  m_element(*m_grid, fem::Q1Quadrature4()),
51  m_rebuild_J_state(true)
52 {
53 
54  PetscErrorCode ierr;
55 
57  m_velocity_shared->metadata(0) = m_velocity.metadata(0);
58  m_velocity_shared->metadata(1) = m_velocity.metadata(1);
59 
60  ierr = DMSetMatType(*m_da, MATBAIJ);
61  PISM_CHK(ierr, "DMSetMatType");
62 
63  ierr = DMCreateMatrix(*m_da, m_J_state.rawptr());
64  PISM_CHK(ierr, "DMCreateMatrix");
65 
66  ierr = KSPCreate(m_grid->com, m_ksp.rawptr());
67  PISM_CHK(ierr, "KSPCreate");
68 
69  double ksp_rtol = 1e-12;
70  ierr = KSPSetTolerances(m_ksp, ksp_rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
71  PISM_CHK(ierr, "KSPSetTolerances");
72 
73  PC pc;
74  ierr = KSPGetPC(m_ksp, &pc);
75  PISM_CHK(ierr, "KSPGetPC");
76 
77  ierr = PCSetType(pc, PCBJACOBI);
78  PISM_CHK(ierr, "PCSetType");
79 
80  ierr = KSPSetFromOptions(m_ksp);
81  PISM_CHK(ierr, "KSPSetFromOptions");
82 }
83 
85 
86  SSAFEM::init();
87 
88  // Get most of the inputs from IceGrid::variables() and fake the rest.
89  //
90  // I will need to fix this at some point.
91  {
92  Geometry geometry(m_grid);
93  geometry.ice_thickness.copy_from(*m_grid->variables().get_2d_scalar("land_ice_thickness"));
94  geometry.bed_elevation.copy_from(*m_grid->variables().get_2d_scalar("bedrock_altitude"));
95  geometry.sea_level_elevation.set(0.0); // FIXME: this should be an input
96  geometry.ice_area_specific_volume.set(0.0);
97 
98  geometry.ensure_consistency(m_config->get_number("stress_balance.ice_free_thickness_standard"));
99 
100  stressbalance::Inputs inputs;
101 
102  inputs.geometry = &geometry;
103  inputs.basal_melt_rate = NULL;
104  inputs.basal_yield_stress = m_grid->variables().get_2d_scalar("tauc");
105  inputs.enthalpy = m_grid->variables().get_3d_scalar("enthalpy");
106  inputs.age = NULL;
107  inputs.bc_mask = m_grid->variables().get_2d_mask("vel_bc_mask");
108  inputs.bc_values = m_grid->variables().get_2d_vector("vel_bc");
109 
110  inputs.water_column_pressure = NULL;
111 
112  cache_inputs(inputs);
113  }
114 }
115 
116 //! Sets the current value of of the design paramter \f$\zeta\f$.
117 /*! This method sets \f$\zeta\f$ but does not solve the %SSA.
118 It it intended for inverse methods that simultaneously compute
119 the pair \f$u\f$ and \f$\zeta\f$ without ever solving the %SSA
120 directly. Use this method in conjuction with
121 \ref assemble_jacobian_state and \ref apply_jacobian_design and their friends.
122 The vector \f$\zeta\f$ is not copied; a reference to the IceModelVec is
123 kept.
124 */
126 
127  m_zeta = &new_zeta;
128 
129  // Convert zeta to hardav.
131 
132  // Cache hardav at the quadrature points.
134 
135  for (PointsWithGhosts p(*m_grid); p; p.next()) {
136  const int i = p.i(), j = p.j();
137  m_coefficients(i, j).hardness = m_hardav(i, j);
138  }
139 
140  // Flag the state jacobian as needing rebuilding.
141  m_rebuild_J_state = true;
142 }
143 
144 //! Sets the current value of the design variable \f$\zeta\f$ and solves the %SSA to find the associated \f$u_{\rm SSA}\f$.
145 /* Use this method for inverse methods employing the reduced gradient. Use this method
146 in conjuction with apply_linearization and apply_linearization_transpose.*/
148  this->set_design(zeta);
149  return this->solve_nocache();
150 }
151 
152 //! Computes the residual function \f$\mathcal{R}(u, \zeta)\f$ as defined in the class-level documentation.
153 /* The value of \f$\zeta\f$ is set prior to this call via set_design or linearize_at. The value
154 of the residual is returned in \a RHS.*/
156 
157  IceModelVec::AccessList l{&u, &RHS};
158 
159  this->compute_local_function(u.array(), RHS.array());
160 }
161 
162 //! Computes the residual function \f$\mathcal{R}(u, \zeta)\f$ defined in the class-level documentation.
163 /* The return value is specified via a Vec for the benefit of certain TAO routines. Otherwise,
164 the method is identical to the assemble_residual returning values as a StateVec (an IceModelVec2V).*/
166 
168 
169  petsc::DMDAVecArray rhs_a(m_da, RHS);
170 
171  this->compute_local_function(u.array(), (Vector2**)rhs_a.get());
172 }
173 
174 //! Assembles the state Jacobian matrix.
175 /* The matrix depends on the current value of the design variable \f$\zeta\f$ and the current
176 value of the state variable \f$u\f$. The specification of \f$\zeta\f$ is done earlier
177 with set_design or linearize_at. The value of \f$u\f$ is specified explicitly as an argument
178 to this method.
179  @param[in] u Current state variable value.
180  @param[out] J computed state Jacobian.
181 */
184 
185  this->compute_local_jacobian(u.array(), Jac);
186 }
187 
188 //! Applies the design Jacobian matrix to a perturbation of the design variable.
189 /*! The return value uses a DesignVector (IceModelVec2V), which can be ghostless. Ghosts (if present) are updated.
190 \overload
191 */
193  IceModelVec2S &dzeta,
194  IceModelVec2V &du) {
196 
197  this->apply_jacobian_design(u, dzeta, du.array());
198 }
199 
200 //! Applies the design Jacobian matrix to a perturbation of the design variable.
201 /*! The return value is a Vec for the benefit of TAO. It is assumed to be ghostless; no communication is done.
202 \overload
203 */
205  IceModelVec2S &dzeta,
206  Vec du) {
207  petsc::DMDAVecArray du_a(m_da, du);
208  this->apply_jacobian_design(u, dzeta, (Vector2**)du_a.get());
209 }
210 
211 //! @brief Applies the design Jacobian matrix to a perturbation of the
212 //! design variable.
213 
214 /*! The matrix depends on the current value of the design variable
215  \f$\zeta\f$ and the current value of the state variable \f$u\f$.
216  The specification of \f$\zeta\f$ is done earlier with set_design
217  or linearize_at. The value of \f$u\f$ is specified explicitly as
218  an argument to this method.
219 
220  @param[in] u Current state variable value.
221 
222  @param[in] dzeta Perturbation of the design variable. Prefers
223  vectors with ghosts; will copy to a ghosted vector
224  if needed.
225 
226  @param[out] du_a Computed corresponding perturbation of the state
227  variable. The array \a du_a should be extracted
228  first from a Vec or an IceModelVec.
229 
230  Typically this method is called via one of its overloads.
231 */
233  IceModelVec2S &dzeta,
234  Vector2 **du_a) {
235 
236  const unsigned int Nk = fem::q1::n_chi;
237  const unsigned int Nq = m_element.n_pts();
238  const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
239 
241 
242  IceModelVec2S *dzeta_local;
243  if (dzeta.stencil_width() > 0) {
244  dzeta_local = &dzeta;
245  } else {
246  m_dzeta_local.copy_from(dzeta);
247  dzeta_local = &m_dzeta_local;
248  }
249  list.add(*dzeta_local);
250 
251  // Zero out the portion of the function we are responsible for computing.
252  for (Points p(*m_grid); p; p.next()) {
253  const int i = p.i(), j = p.j();
254 
255  du_a[j][i].u = 0.0;
256  du_a[j][i].v = 0.0;
257  }
258 
259  // Aliases to help with notation consistency below.
260  const IceModelVec2Int *dirichletLocations = &m_bc_mask;
261  const IceModelVec2V *dirichletValues = &m_bc_values;
262  double dirichletWeight = m_dirichletScale;
263 
264  Vector2 u_e[Nk];
265  Vector2 U[Nq_max], U_x[Nq_max], U_y[Nq_max];
266 
267  Vector2 du_e[Nk];
268 
269  double dzeta_e[Nk];
270 
271  double zeta_e[Nk];
272 
273  double dB_e[Nk];
274  double dB_q[Nq_max];
275 
276  fem::DirichletData_Vector dirichletBC(dirichletLocations, dirichletValues,
277  dirichletWeight);
279 
280  // Loop through all elements.
281  const int
282  xs = m_element_index.xs,
283  xm = m_element_index.xm,
284  ys = m_element_index.ys,
285  ym = m_element_index.ym;
286 
287  ParallelSection loop(m_grid->com);
288  try {
289  for (int j =ys; j<ys+ym; j++) {
290  for (int i =xs; i<xs+xm; i++) {
291 
292  // Zero out the element-local residual in prep for updating it.
293  for (unsigned int k=0; k<Nk; k++) {
294  du_e[k].u = 0;
295  du_e[k].v = 0;
296  }
297 
298  // Initialize the map from global to local degrees of freedom for this element.
299  m_element.reset(i, j);
300 
301  // Obtain the value of the solution at the nodes adjacent to the element,
302  // fix dirichlet values, and compute values at quad pts.
303  m_element.nodal_values(u.array(), u_e);
304  if (dirichletBC) {
305  dirichletBC.constrain(m_element);
306  dirichletBC.enforce(m_element, u_e);
307  }
308  m_element.evaluate(u_e, U, U_x, U_y);
309 
310  // Compute dzeta at the nodes
311  m_element.nodal_values(dzeta_local->array(), dzeta_e);
312  if (fixedZeta) {
313  fixedZeta.enforce_homogeneous(m_element, dzeta_e);
314  }
315 
316  // Compute the change in hardav with respect to zeta at the quad points.
317  m_element.nodal_values(m_zeta->array(), zeta_e);
318  for (unsigned int k=0; k<Nk; k++) {
319  m_design_param.toDesignVariable(zeta_e[k], NULL, dB_e + k);
320  dB_e[k]*=dzeta_e[k];
321  }
322  m_element.evaluate(dB_e, dB_q);
323 
324  double thickness[Nq_max];
325  {
326  Coefficients coeffs[Nk];
327  int mask[Nq_max];
328  double tauc[Nq_max];
329  double hardness[Nq_max];
330 
331  m_element.nodal_values(m_coefficients.array(), coeffs);
332 
334  mask, thickness, tauc, hardness);
335  }
336 
337  for (unsigned int q = 0; q < Nq; q++) {
338  // Symmetric gradient at the quadrature point.
339  double Duqq[3] = {U_x[q].u, U_y[q].v, 0.5 * (U_y[q].u + U_x[q].v)};
340 
341  double d_nuH = 0;
342  if (thickness[q] >= strength_extension->get_min_thickness()) {
343  m_flow_law->effective_viscosity(dB_q[q],
344  secondInvariant_2D(U_x[q], U_y[q]),
345  &d_nuH, NULL);
346  d_nuH *= (2.0 * thickness[q]);
347  }
348 
349  auto W = m_element.weight(q);
350 
351  for (unsigned int k = 0; k < Nk; k++) {
352  const fem::Germ &testqk = m_element.chi(q, k);
353  du_e[k].u += W*d_nuH*(testqk.dx*(2*Duqq[0] + Duqq[1]) + testqk.dy*Duqq[2]);
354  du_e[k].v += W*d_nuH*(testqk.dy*(2*Duqq[1] + Duqq[0]) + testqk.dx*Duqq[2]);
355  }
356  } // q
357  m_element.add_contribution(du_e, du_a);
358  } // j
359  } // i
360  } catch (...) {
361  loop.failed();
362  }
363  loop.check();
364 
365  if (dirichletBC) {
366  dirichletBC.fix_residual_homogeneous(du_a);
367  }
368 }
369 
370 //! Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
371 /*! The return value uses a StateVector (IceModelVec2S) which can be ghostless; ghosts (if present) are updated.
372 \overload
373 */
375  IceModelVec2V &du,
376  IceModelVec2S &dzeta) {
377  IceModelVec::AccessList l{&dzeta};
378  this->apply_jacobian_design_transpose(u, du, dzeta.array());
379 }
380 
381 //! Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
382 /*! The return value uses a Vec for the benefit of TAO. It is assumed to be ghostless; no communication is done.
383 \overload */
385  IceModelVec2V &du,
386  Vec dzeta) {
387 
388  petsc::DM::Ptr da2 = m_grid->get_dm(1, m_config->get_number("grid.max_stencil_width"));
389  petsc::DMDAVecArray dzeta_a(da2, dzeta);
390  this->apply_jacobian_design_transpose(u, du, (double**)dzeta_a.get());
391 }
392 
393 //! @brief Applies the transpose of the design Jacobian matrix to a
394 //! perturbation of the state variable.
395 
396 /*! The matrix depends on the current value of the design variable
397  \f$\zeta\f$ and the current value of the state variable \f$u\f$.
398  The specification of \f$\zeta\f$ is done earlier with set_design
399  or linearize_at. The value of \f$u\f$ is specified explicitly as
400  an argument to this method.
401 
402  @param[in] u Current state variable value.
403 
404  @param[in] du Perturbation of the state variable. Prefers vectors
405  with ghosts; will copy to a ghosted vector if need be.
406 
407  @param[out] dzeta_a Computed corresponding perturbation of the
408  design variable. The array \a dzeta_a should be
409  extracted first from a Vec or an IceModelVec.
410 
411  Typically this method is called via one of its overloads.
412 */
414  IceModelVec2V &du,
415  double **dzeta_a) {
416 
417  const unsigned int Nk = fem::q1::n_chi;
418  const unsigned int Nq = m_element.n_pts();
419  const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
420 
422 
423  IceModelVec2V *du_local;
424  if (du.stencil_width() > 0) {
425  du_local = &du;
426  } else {
427  m_du_local.copy_from(du);
428  du_local = &m_du_local;
429  }
430  list.add(*du_local);
431 
432  Vector2 u_e[Nk];
433  Vector2 U[Nq_max], U_x[Nq_max], U_y[Nq_max];
434 
435  Vector2 du_e[Nk];
436  Vector2 du_q[Nq_max];
437  Vector2 du_dx_q[Nq_max];
438  Vector2 du_dy_q[Nq_max];
439 
440  double dzeta_e[Nk];
441 
442  // Aliases to help with notation consistency.
443  const IceModelVec2Int *dirichletLocations = &m_bc_mask;
444  const IceModelVec2V *dirichletValues = &m_bc_values;
445  double dirichletWeight = m_dirichletScale;
446 
447  fem::DirichletData_Vector dirichletBC(dirichletLocations, dirichletValues,
448  dirichletWeight);
449 
450  // Zero out the portion of the function we are responsible for computing.
451  for (Points p(*m_grid); p; p.next()) {
452  const int i = p.i(), j = p.j();
453 
454  dzeta_a[j][i] = 0;
455  }
456 
457  const int
458  xs = m_element_index.xs,
459  xm = m_element_index.xm,
460  ys = m_element_index.ys,
461  ym = m_element_index.ym;
462 
463  ParallelSection loop(m_grid->com);
464  try {
465  for (int j = ys; j < ys + ym; j++) {
466  for (int i = xs; i < xs + xm; i++) {
467  // Initialize the map from global to local degrees of freedom for this element.
468  m_element.reset(i, j);
469 
470  // Obtain the value of the solution at the nodes adjacent to the element.
471  // Compute the solution values and symmetric gradient at the quadrature points.
472  m_element.nodal_values(du.array(), du_e);
473  if (dirichletBC) {
474  dirichletBC.enforce_homogeneous(m_element, du_e);
475  }
476  m_element.evaluate(du_e, du_q, du_dx_q, du_dy_q);
477 
478  m_element.nodal_values(u.array(), u_e);
479  if (dirichletBC) {
480  dirichletBC.enforce(m_element, u_e);
481  }
482  m_element.evaluate(u_e, U, U_x, U_y);
483 
484  // Zero out the element-local residual in prep for updating it.
485  for (unsigned int k = 0; k < Nk; k++) {
486  dzeta_e[k] = 0;
487  }
488 
489  double thickness[Nq_max];
490  {
491  Coefficients coeffs[Nk];
492  int mask[Nq_max];
493  double tauc[Nq_max];
494  double hardness[Nq_max];
495 
496  m_element.nodal_values(m_coefficients.array(), coeffs);
497 
499  mask, thickness, tauc, hardness);
500  }
501 
502  for (unsigned int q = 0; q < Nq; q++) {
503  // Symmetric gradient at the quadrature point.
504  double Duqq[3] = {U_x[q].u, U_y[q].v, 0.5 * (U_y[q].u + U_x[q].v)};
505 
506  // Determine "d_nuH / dB" at the quadrature point
507  double d_nuH_dB = 0;
508  if (thickness[q] >= strength_extension->get_min_thickness()) {
509  m_flow_law->effective_viscosity(1.0,
510  secondInvariant_2D(U_x[q], U_y[q]),
511  &d_nuH_dB, NULL);
512  d_nuH_dB *= (2.0 * thickness[q]);
513  }
514 
515  auto W = m_element.weight(q);
516 
517  for (unsigned int k = 0; k < Nk; k++) {
518  dzeta_e[k] += W*d_nuH_dB*m_element.chi(q, k).val*((du_dx_q[q].u*(2*Duqq[0] + Duqq[1]) +
519  du_dy_q[q].u*Duqq[2]) +
520  (du_dy_q[q].v*(2*Duqq[1] + Duqq[0]) +
521  du_dx_q[q].v*Duqq[2]));
522  }
523  } // q
524 
525  m_element.add_contribution(dzeta_e, dzeta_a);
526  } // j
527  } // i
528  } catch (...) {
529  loop.failed();
530  }
531  loop.check();
532 
533  for (Points p(*m_grid); p; p.next()) {
534  const int i = p.i(), j = p.j();
535 
536  double dB_dzeta;
537  m_design_param.toDesignVariable((*m_zeta)(i, j), NULL, &dB_dzeta);
538  dzeta_a[j][i] *= dB_dzeta;
539  }
540 
543  fixedZeta.fix_residual_homogeneous(dzeta_a);
544  }
545 }
546 
547 /*!\brief Applies the linearization of the forward map (i.e. the reduced gradient \f$DF\f$ described in
548 the class-level documentation.) */
549 /*! As described previously,
550 \f[
551 Df = J_{\rm State}^{-1} J_{\rm Design}.
552 \f]
553 Applying the linearization then involves the solution of a linear equation.
554 The matrices \f$J_{\rm State}\f$ and \f$J_{\rm Design}\f$ both depend on the value of the
555 design variable \f$\zeta\f$ and the value of the corresponding state variable \f$u=F(\zeta)\f$.
556 These are established by first calling linearize_at.
557  @param[in] dzeta Perturbation of the design variable
558  @param[out] du Computed corresponding perturbation of the state variable; ghosts (if present) are updated.
559 */
561 
562  PetscErrorCode ierr;
563 
564  if (m_rebuild_J_state) {
566  m_rebuild_J_state = false;
567  }
568 
570  m_du_global.scale(-1);
571 
572  // call PETSc to solve linear system by iterative method.
573  ierr = KSPSetOperators(m_ksp, m_J_state, m_J_state);
574  PISM_CHK(ierr, "KSPSetOperators");
575 
576  ierr = KSPSolve(m_ksp, m_du_global.vec(), m_du_global.vec());
577  PISM_CHK(ierr, "KSPSolve"); // SOLVE
578 
579  KSPConvergedReason reason;
580  ierr = KSPGetConvergedReason(m_ksp, &reason);
581  PISM_CHK(ierr, "KSPGetConvergedReason");
582 
583  if (reason < 0) {
584  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IP_SSAHardavForwardProblem::apply_linearization solve"
585  " failed to converge (KSP reason %s)",
586  KSPConvergedReasons[reason]);
587  } else {
588  m_log->message(4,
589  "IP_SSAHardavForwardProblem::apply_linearization converged"
590  " (KSP reason %s)\n",
591  KSPConvergedReasons[reason]);
592  }
593 
595 }
596 
597 /*! \brief Applies the transpose of the linearization of the forward map
598  (i.e. the transpose of the reduced gradient \f$DF\f$ described in the class-level documentation.) */
599 /*! As described previously,
600 \f[
601 Df = J_{\rm State}^{-1} J_{\rm Design}.
602 \f]
603 so
604 \f[
605 Df^t = J_{\rm Design}^t \; (J_{\rm State}^t)^{-1} .
606 \f]
607 Applying the transpose of the linearization then involves the solution of a linear equation.
608 The matrices \f$J_{\rm State}\f$ and \f$J_{\rm Design}\f$ both depend on the value of the
609 design variable \f$\zeta\f$ and the value of the corresponding state variable \f$u=F(\zeta)\f$.
610 These are established by first calling linearize_at.
611  @param[in] du Perturbation of the state variable
612  @param[out] dzeta Computed corresponding perturbation of the design variable; ghosts (if present) are updated.
613 */
615  IceModelVec2S &dzeta) {
616 
617  PetscErrorCode ierr;
618 
619  if (m_rebuild_J_state) {
621  m_rebuild_J_state = false;
622  }
623 
624  // Aliases to help with notation consistency below.
625  const IceModelVec2Int *dirichletLocations = &m_bc_mask;
626  const IceModelVec2V *dirichletValues = &m_bc_values;
627  double dirichletWeight = m_dirichletScale;
628 
630 
631  fem::DirichletData_Vector dirichletBC(dirichletLocations, dirichletValues,
632  dirichletWeight);
633  if (dirichletBC) {
636  }
637 
638  // call PETSc to solve linear system by iterative method.
639  ierr = KSPSetOperators(m_ksp, m_J_state, m_J_state);
640  PISM_CHK(ierr, "KSPSetOperators");
641 
642  ierr = KSPSolve(m_ksp, m_du_global.vec(), m_du_global.vec());
643  PISM_CHK(ierr, "KSPSolve"); // SOLVE
644 
645  KSPConvergedReason reason;
646  ierr = KSPGetConvergedReason(m_ksp, &reason);
647  PISM_CHK(ierr, "KSPGetConvergedReason");
648 
649  if (reason < 0) {
650  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IP_SSAHardavForwardProblem::apply_linearization solve failed to converge (KSP reason %s)",
651  KSPConvergedReasons[reason]);
652  } else {
653  m_log->message(4,
654  "IP_SSAHardavForwardProblem::apply_linearization converged (KSP reason %s)\n",
655  KSPConvergedReasons[reason]);
656  }
657 
659  dzeta.scale(-1);
660 
661  if (dzeta.stencil_width() > 0) {
662  dzeta.update_ghosts();
663  }
664 }
665 
666 } // end of namespace inverse
667 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
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
void ensure_consistency(double ice_free_thickness_threshold)
Definition: Geometry.cc:122
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S ice_area_specific_volume
Definition: Geometry.hh:54
IceModelVec2S sea_level_elevation
Definition: Geometry.hh:50
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
void add(double alpha, const IceModelVec2S &x)
void copy_from(const IceModelVec2S &source)
void copy_from(const IceModelVec2< T > &source)
Definition: IceModelVec2.hh:97
void add(double alpha, const IceModelVec2< T > &x)
Definition: IceModelVec2.hh:89
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:334
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
petsc::Vec & vec() const
Definition: iceModelVec.cc:342
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
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_homogeneous(const Element2 &element, double *x_e)
void fix_residual_homogeneous(double **r_global)
void enforce(const Element2 &element, Vector2 *x_e)
void enforce_homogeneous(const Element2 &element, Vector2 *x_e)
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
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
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
double weight(unsigned int q) const
Weight of the quadrature point q
Definition: Element.hh:85
virtual void convertToDesignVariable(IceModelVec2S &zeta, IceModelVec2S &d, bool communicate=true)
Transforms a vector of values to a vector of values.
virtual void toDesignVariable(double zeta, double *value, double *derivative)=0
Converts from parameterization value to .
petsc::Mat m_J_state
Mat used in apply_linearization and apply_linearization_transpose.
virtual void apply_linearization_transpose(IceModelVec2V &du, IceModelVec2S &dzeta)
Applies the transpose of the linearization of the forward map (i.e. the transpose of the reduced grad...
IPDesignVariableParameterization & m_design_param
The function taking to .
virtual TerminationReason::Ptr linearize_at(IceModelVec2S &zeta)
Sets the current value of the design variable and solves the SSA to find the associated .
IceModelVec2S * m_zeta
Current value of zeta, provided from caller.
virtual void set_design(IceModelVec2S &zeta)
Sets the current value of of the design paramter .
virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, IceModelVec2V &du)
Applies the design Jacobian matrix to a perturbation of the design variable.
virtual void assemble_residual(IceModelVec2V &u, IceModelVec2V &R)
Computes the residual function as defined in the class-level documentation.
IceModelVec2Int * m_fixed_design_locations
Locations where should not be adjusted.
IceModelVec2V m_du_global
Temporary storage when state vectors need to be used without ghosts.
virtual void assemble_jacobian_state(IceModelVec2V &u, Mat J)
Assembles the state Jacobian matrix.
bool m_rebuild_J_state
Flag indicating that the state jacobian matrix needs rebuilding.
IP_SSAHardavForwardProblem(IceGrid::ConstPtr g, IPDesignVariableParameterization &tp)
Constructs from the same objects as SSAFEM, plus a specification of how is parameterized.
IceModelVec2S m_dzeta_local
Storage for d_zeta with ghosts, if needed when an argument d_zeta is ghost-less.
virtual void apply_linearization(IceModelVec2S &dzeta, IceModelVec2V &du)
Applies the linearization of the forward map (i.e. the reduced gradient described in the class-level...
IceModelVec2V m_du_local
Temporary storage when state vectors need to be used with ghosts.
IceModelVec2S m_hardav
Vertically-averaged ice hardness.
virtual void apply_jacobian_design_transpose(IceModelVec2V &u, IceModelVec2V &du, IceModelVec2S &dzeta)
Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
petsc::KSP m_ksp
KSP used in apply_linearization and apply_linearization_transpose.
std::shared_ptr< Wrapper > Ptr
Definition: Wrapper.hh:31
const IceModelVec2S * water_column_pressure
const IceModelVec3 * enthalpy
const IceModelVec2S * basal_yield_stress
const IceModelVec3 * age
const IceModelVec2V * bc_values
const IceModelVec2Int * bc_mask
const IceModelVec2S * basal_melt_rate
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
void cache_inputs(const Inputs &inputs)
Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve.
Definition: SSAFEM.cc:261
IceModelVec2< Coefficients > m_coefficients
Definition: SSAFEM.hh:76
void compute_local_function(Vector2 const *const *const velocity, Vector2 **residual)
Implements the callback for computing the residual.
Definition: SSAFEM.cc:729
IceModelVec2Int m_bc_mask
Definition: SSAFEM.hh:69
TerminationReason::Ptr solve_nocache()
Definition: SSAFEM.cc:187
void compute_local_jacobian(Vector2 const *const *const velocity, Mat J)
Implements the callback for computing the Jacobian.
Definition: SSAFEM.cc:951
IceModelVec2V m_bc_values
Definition: SSAFEM.hh:70
double get_min_thickness() const
Returns minimum thickness to trigger use of extension.
Definition: SSA.cc:70
SSAStrengthExtension * strength_extension
Definition: SSA.hh:115
std::shared_ptr< petsc::DM > m_da
Definition: SSA.hh:146
std::shared_ptr< rheology::FlowLaw > m_flow_law
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
const int n_chi
Definition: FEM.hh:191
const unsigned int MAX_QUADRATURE_SIZE
Definition: FEM.hh:173
static const double g
Definition: exactTestP.cc:39
static const double k
Definition: exactTestP.cc:45
static double secondInvariant_2D(const Vector2 &U_x, const Vector2 &U_y)
Definition: FlowLaw.hh:44
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
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