tIP_SSATaucForwardProblem.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tIP_SSATaucForwardProblem.cc (24047B)
---
1 // Copyright (C) 2012, 2014, 2015, 2016, 2017, 2019 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
19 #include "IP_SSATaucForwardProblem.hh"
20 #include "pism/basalstrength/basal_resistance.hh"
21 #include "pism/util/IceGrid.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/pism_utilities.hh"
26 #include "pism/geometry/Geometry.hh"
27 #include "pism/stressbalance/StressBalance.hh"
28
29 namespace pism {
30 namespace inverse {
31
32 IP_SSATaucForwardProblem::IP_SSATaucForwardProblem(IceGrid::ConstPtr g,
33 IPDesignVariableParameterization &tp)
34 : SSAFEM(g),
35 m_zeta(NULL),
36 m_fixed_tauc_locations(NULL),
37 m_tauc_param(tp),
38 m_element_index(*m_grid),
39 m_element(*m_grid),
40 m_quadrature(g->dx(), g->dy(), 1.0),
41 m_rebuild_J_state(true) {
42
43 PetscErrorCode ierr;
44 int stencil_width = 1;
45
46 m_velocity_shared.reset(new IceModelVec2V(m_grid, "dummy", WITHOUT_GHOSTS));
47 m_velocity_shared->metadata(0) = m_velocity.metadata(0);
48 m_velocity_shared->metadata(1) = m_velocity.metadata(1);
49
50 m_dzeta_local.create(m_grid, "d_zeta_local", WITH_GHOSTS, stencil_width);
51
52 m_du_global.create(m_grid, "linearization work vector (sans ghosts)", WITHOUT_GHOSTS, stencil_width);
53 m_du_local.create(m_grid, "linearization work vector (with ghosts)", WITH_GHOSTS, stencil_width);
54
55 m_tauc_copy.create(m_grid, "tauc", WITH_GHOSTS, m_config->get_number("grid.max_stencil_width"));
56 m_tauc_copy.set_attrs("diagnostic",
57 "yield stress for basal till (plastic or pseudo-plastic model)",
58 "Pa", "Pa", "", 0);
59
60 ierr = DMSetMatType(*m_da, MATBAIJ);
61 PISM_CHK(ierr, "DMSetMatType");
62 ierr = DMCreateMatrix(*m_da, m_J_state.rawptr());
63 PISM_CHK(ierr, "DMCreateMatrix");
64
65 ierr = KSPCreate(m_grid->com, m_ksp.rawptr());
66 PISM_CHK(ierr, "KSPCreate");
67
68 double ksp_rtol = 1e-12;
69 ierr = KSPSetTolerances(m_ksp, ksp_rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
70 PISM_CHK(ierr, "KSPSetTolerances");
71
72 PC pc;
73 ierr = KSPGetPC(m_ksp, &pc);
74 PISM_CHK(ierr, "KSPGetPC");
75
76 ierr = PCSetType(pc, PCBJACOBI);
77 PISM_CHK(ierr, "PCSetType");
78
79 ierr = KSPSetFromOptions(m_ksp);
80 PISM_CHK(ierr, "KSPSetFromOptions");
81 }
82
83 IP_SSATaucForwardProblem::~IP_SSATaucForwardProblem() {
84 // empty
85 }
86
87 void IP_SSATaucForwardProblem::init() {
88
89 // This calls SSA::init(), which calls pism::Vars::get_2d_scalar()
90 // to set m_tauc.
91 SSAFEM::init();
92
93 // Get most of the inputs from IceGrid::variables() and fake the rest.
94 //
95 // I will need to fix this at some point.
96 {
97 Geometry geometry(m_grid);
98 geometry.ice_thickness.copy_from(*m_grid->variables().get_2d_scalar("land_ice_thickness"));
99 geometry.bed_elevation.copy_from(*m_grid->variables().get_2d_scalar("bedrock_altitude"));
100 geometry.sea_level_elevation.set(0.0);
101 geometry.ice_area_specific_volume.set(0.0);
102
103 geometry.ensure_consistency(m_config->get_number("stress_balance.ice_free_thickness_standard"));
104
105 stressbalance::Inputs inputs;
106
107 inputs.geometry = &geometry;
108 inputs.basal_melt_rate = NULL;
109 inputs.melange_back_pressure = NULL;
110 inputs.basal_yield_stress = m_grid->variables().get_2d_scalar("tauc");
111 inputs.enthalpy = m_grid->variables().get_3d_scalar("enthalpy");
112 inputs.age = NULL;
113 inputs.bc_mask = m_grid->variables().get_2d_mask("bc_mask");
114 inputs.bc_values = m_grid->variables().get_2d_vector("vel_ssa_bc");
115
116 cache_inputs(inputs);
117 }
118 }
119
120 //! Sets the current value of of the design paramter \f$\zeta\f$.
121 /*! This method sets \f$\zeta\f$ but does not solve the %SSA.
122 It it intended for inverse methods that simultaneously compute
123 the pair \f$u\f$ and \f$\zeta\f$ without ever solving the %SSA
124 directly. Use this method in conjuction with
125 \ref assemble_jacobian_state and \ref apply_jacobian_design and their friends.
126 The vector \f$\zeta\f$ is not copied; a reference to the IceModelVec is
127 kept.
128 */
129 void IP_SSATaucForwardProblem::set_design(IceModelVec2S &new_zeta) {
130
131 IceModelVec2S &tauc = m_tauc_copy;
132
133 m_zeta = &new_zeta;
134
135 // Convert zeta to tauc.
136 m_tauc_param.convertToDesignVariable(*m_zeta, tauc);
137
138 // Cache tauc at the quadrature points.
139 IceModelVec::AccessList list{&tauc, &m_coefficients};
140
141 for (PointsWithGhosts p(*m_grid); p; p.next()) {
142 const int i = p.i(), j = p.j();
143 m_coefficients(i, j).tauc = tauc(i, j);
144 }
145
146 // Flag the state jacobian as needing rebuilding.
147 m_rebuild_J_state = true;
148 }
149
150 //! Sets the current value of the design variable \f$\zeta\f$ and solves the %SSA to find the associated \f$u_{\rm SSA}\f$.
151 /* Use this method for inverse methods employing the reduced gradient. Use this method
152 in conjuction with apply_linearization and apply_linearization_transpose.*/
153 TerminationReason::Ptr IP_SSATaucForwardProblem::linearize_at(IceModelVec2S &zeta) {
154 this->set_design(zeta);
155 return this->solve_nocache();
156 }
157
158 //! Computes the residual function \f$\mathcal{R}(u, \zeta)\f$ as defined in the class-level documentation.
159 /* The value of \f$\zeta\f$ is set prior to this call via set_design or linearize_at. The value
160 of the residual is returned in \a RHS.*/
161 void IP_SSATaucForwardProblem::assemble_residual(IceModelVec2V &u, IceModelVec2V &RHS) {
162
163 Vector2
164 **u_a = u.get_array(),
165 **rhs_a = RHS.get_array();
166
167 this->compute_local_function(u_a, rhs_a);
168
169 u.end_access();
170 RHS.end_access();
171 }
172
173 //! Computes the residual function \f$\mathcal{R}(u, \zeta)\f$ defined in the class-level documentation.
174 /* The return value is specified via a Vec for the benefit of certain TAO routines. Otherwise,
175 the method is identical to the assemble_residual returning values as a StateVec (an IceModelVec2V).*/
176 void IP_SSATaucForwardProblem::assemble_residual(IceModelVec2V &u, Vec RHS) {
177
178 Vector2 **u_a = u.get_array();
179 petsc::DMDAVecArray rhs_a(m_da, RHS);
180 this->compute_local_function(u_a, (Vector2**)rhs_a.get());
181 u.end_access();
182 }
183
184 //! Assembles the state Jacobian matrix.
185 /* The matrix depends on the current value of the design variable \f$\zeta\f$ and the current
186 value of the state variable \f$u\f$. The specification of \f$\zeta\f$ is done earlier
187 with set_design or linearize_at. The value of \f$u\f$ is specified explicitly as an argument
188 to this method.
189 @param[in] u Current state variable value.
190 @param[out] J computed state Jacobian.
191 */
192 void IP_SSATaucForwardProblem::assemble_jacobian_state(IceModelVec2V &u, Mat Jac) {
193
194 Vector2 **u_a = u.get_array();
195
196 this->compute_local_jacobian(u_a, Jac);
197
198 u.end_access();
199 }
200
201 //! Applies the design Jacobian matrix to a perturbation of the design variable.
202 /*! The return value uses a DesignVector (IceModelVec2V), which can be ghostless. Ghosts (if present) are updated.
203 \overload
204 */
205 void IP_SSATaucForwardProblem::apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta,
206 IceModelVec2V &du) {
207 Vector2 **du_a = du.get_array();
208 this->apply_jacobian_design(u, dzeta, du_a);
209 du.end_access();
210 }
211
212 //! Applies the design Jacobian matrix to a perturbation of the design variable.
213 /*! The return value is a Vec for the benefit of TAO. It is assumed to
214 be ghostless; no communication is done. \overload
215 */
216 void IP_SSATaucForwardProblem::apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta,
217 Vec du) {
218 petsc::DMDAVecArray du_a(m_da, du);
219 this->apply_jacobian_design(u, dzeta, (Vector2**)du_a.get());
220 }
221
222 //! Applies the design Jacobian matrix to a perturbation of the design variable.
223 /*! The matrix depends on the current value of the design variable \f$\zeta\f$ and the current
224 value of the state variable \f$u\f$. The specification of \f$\zeta\f$ is done earlier
225 with set_design or linearize_at. The value of \f$u\f$ is specified explicitly as an argument
226 to this method.
227 @param[in] u Current state variable value.
228 @param[in] dzeta Perturbation of the design variable. Prefers vectors with ghosts; will copy to a ghosted vector if needed.
229 @param[out] du_a Computed corresponding perturbation of the state variable. The array \a du_a
230 should be extracted first from a Vec or an IceModelVec.
231
232 Typically this method is called via one of its overloads.
233 */
234 void IP_SSATaucForwardProblem::apply_jacobian_design(IceModelVec2V &u,
235 IceModelVec2S &dzeta,
236 Vector2 **du_a) {
237 const unsigned int Nk = fem::q1::n_chi;
238 const unsigned int Nq = m_quadrature.n();
239 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
240
241 IceModelVec::AccessList list{&m_coefficients, m_zeta, &u};
242
243 IceModelVec2S *dzeta_local;
244 if (dzeta.stencil_width() > 0) {
245 dzeta_local = &dzeta;
246 } else {
247 m_dzeta_local.copy_from(dzeta);
248 dzeta_local = &m_dzeta_local;
249 }
250 list.add(*dzeta_local);
251
252 // Zero out the portion of the function we are responsible for computing.
253 for (Points p(*m_grid); p; p.next()) {
254 const int i = p.i(), j = p.j();
255
256 du_a[j][i].u = 0.0;
257 du_a[j][i].v = 0.0;
258 }
259
260 // Aliases to help with notation consistency below.
261 const IceModelVec2Int *dirichletLocations = m_bc_mask;
262 const IceModelVec2V *dirichletValues = m_bc_values;
263 double dirichletWeight = m_dirichletScale;
264
265 Vector2 u_e[Nk];
266 Vector2 u_q[Nq_max];
267
268 Vector2 du_e[Nk];
269
270 double dzeta_e[Nk];
271
272 double zeta_e[Nk];
273
274 double dtauc_e[Nk];
275 double dtauc_q[Nq_max];
276
277 // An Nq by Nk array of test function values.
278 const fem::Germs *test = m_quadrature.test_function_values();
279
280 fem::DirichletData_Vector dirichletBC(dirichletLocations, dirichletValues,
281 dirichletWeight);
282 fem::DirichletData_Scalar fixedZeta(m_fixed_tauc_locations, NULL);
283
284 // Jacobian times weights for quadrature.
285 const double* W = m_quadrature.weights();
286
287 // Loop through all elements.
288 const int
289 xs = m_element_index.xs,
290 xm = m_element_index.xm,
291 ys = m_element_index.ys,
292 ym = m_element_index.ym;
293
294 ParallelSection loop(m_grid->com);
295 try {
296 for (int j = ys; j < ys + ym; j++) {
297 for (int i = xs; i < xs + xm; i++) {
298
299 // Zero out the element residual in prep for updating it.
300 for (unsigned int k = 0; k < Nk; k++) {
301 du_e[k].u = 0;
302 du_e[k].v = 0;
303 }
304
305 // Initialize the map from global to local degrees of freedom for this element.
306 m_element.reset(i, j);
307
308 // Obtain the value of the solution at the nodes adjacent to the element,
309 // fix dirichlet values, and compute values at quad pts.
310 m_element.nodal_values(u, u_e);
311 if (dirichletBC) {
312 dirichletBC.constrain(m_element);
313 dirichletBC.enforce(m_element, u_e);
314 }
315 quadrature_point_values(m_quadrature, u_e, u_q);
316
317 // Compute dzeta at the nodes
318 m_element.nodal_values(*dzeta_local, dzeta_e);
319 if (fixedZeta) {
320 fixedZeta.enforce_homogeneous(m_element, dzeta_e);
321 }
322
323 // Compute the change in tau_c with respect to zeta at the quad points.
324 m_element.nodal_values(*m_zeta, zeta_e);
325 for (unsigned int k = 0; k < Nk; k++) {
326 m_tauc_param.toDesignVariable(zeta_e[k], NULL, dtauc_e + k);
327 dtauc_e[k] *= dzeta_e[k];
328 }
329 quadrature_point_values(m_quadrature, dtauc_e, dtauc_q);
330
331 int mask[Nq_max];
332 {
333 Coefficients coeffs[Nk];
334 double thickness[Nq_max];
335 double tauc[Nq_max];
336 double hardness[Nq_max];
337
338 m_element.nodal_values(m_coefficients, coeffs);
339
340 quad_point_values(m_quadrature, coeffs,
341 mask, thickness, tauc, hardness);
342 }
343
344 for (unsigned int q = 0; q < Nq; q++) {
345 Vector2 u_qq = u_q[q];
346
347 // Determine "dbeta / dzeta" at the quadrature point
348 double dbeta = 0;
349 if (mask::grounded_ice(mask[q])) {
350 dbeta = m_basal_sliding_law->drag(dtauc_q[q], u_qq.u, u_qq.v);
351 }
352
353 for (unsigned int k = 0; k < Nk; k++) {
354 du_e[k].u += W[q]*dbeta*u_qq.u*test[q][k].val;
355 du_e[k].v += W[q]*dbeta*u_qq.v*test[q][k].val;
356 }
357 } // q
358 m_element.add_contribution(du_e, du_a);
359 } // j
360 } // i
361 } catch (...) {
362 loop.failed();
363 }
364 loop.check();
365
366 if (dirichletBC) {
367 dirichletBC.fix_residual_homogeneous(du_a);
368 }
369 }
370
371 //! Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
372 /*! The return value uses a StateVector (IceModelVec2S) which can be ghostless; ghosts (if present) are updated.
373 \overload
374 */
375 void IP_SSATaucForwardProblem::apply_jacobian_design_transpose(IceModelVec2V &u,
376 IceModelVec2V &du,
377 IceModelVec2S &dzeta) {
378 double **dzeta_a = dzeta.get_array();
379 this->apply_jacobian_design_transpose(u, du, dzeta_a);
380 dzeta.end_access();
381 }
382
383 //! Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
384 /*! The return value uses a Vec for the benefit of TAO. It is assumed to be ghostless; no communication is done.
385 \overload */
386 void IP_SSATaucForwardProblem::apply_jacobian_design_transpose(IceModelVec2V &u,
387 IceModelVec2V &du,
388 Vec dzeta) {
389 petsc::DM::Ptr da2 = m_grid->get_dm(1, m_config->get_number("grid.max_stencil_width"));
390 petsc::DMDAVecArray dzeta_a(da2, dzeta);
391 this->apply_jacobian_design_transpose(u, du, (double**)dzeta_a.get());
392 }
393
394 //! Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
395 /*! The matrix depends on the current value of the design variable \f$\zeta\f$ and the current
396 value of the state variable \f$u\f$. The specification of \f$\zeta\f$ is done earlier
397 with set_design or linearize_at. The value of \f$u\f$ is specified explicitly as an argument
398 to this method.
399 @param[in] u Current state variable value.
400 @param[in] du Perturbation of the state variable. Prefers vectors with ghosts; will copy to a ghosted vector if need be.
401 @param[out] dzeta_a Computed corresponding perturbation of the design variable. The array \a dzeta_a
402 should be extracted first from a Vec or an IceModelVec.
403
404 Typically this method is called via one of its overloads.
405 */
406 void IP_SSATaucForwardProblem::apply_jacobian_design_transpose(IceModelVec2V &u,
407 IceModelVec2V &du,
408 double **dzeta_a) {
409 const unsigned int Nk = fem::q1::n_chi;
410 const unsigned int Nq = m_quadrature.n();
411 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
412
413 IceModelVec::AccessList list{&m_coefficients, m_zeta, &u};
414
415 IceModelVec2V *du_local;
416 if (du.stencil_width() > 0) {
417 du_local = &du;
418 } else {
419 m_du_local.copy_from(du);
420 du_local = &m_du_local;
421 }
422 list.add(*du_local);
423
424 Vector2 u_e[Nk];
425 Vector2 u_q[Nq_max];
426
427 Vector2 du_e[Nk];
428 Vector2 du_q[Nq_max];
429
430 double dzeta_e[Nk];
431
432 // An Nq by Nk array of test function values.
433 const fem::Germs *test = m_quadrature.test_function_values();
434
435 // Aliases to help with notation consistency.
436 const IceModelVec2Int *dirichletLocations = m_bc_mask;
437 const IceModelVec2V *dirichletValues = m_bc_values;
438 double dirichletWeight = m_dirichletScale;
439
440 fem::DirichletData_Vector dirichletBC(dirichletLocations, dirichletValues,
441 dirichletWeight);
442
443 // Jacobian times weights for quadrature.
444 const double* W = m_quadrature.weights();
445
446 // Zero out the portion of the function we are responsible for computing.
447 for (Points p(*m_grid); p; p.next()) {
448 const int i = p.i(), j = p.j();
449
450 dzeta_a[j][i] = 0;
451 }
452
453 const int
454 xs = m_element_index.xs,
455 xm = m_element_index.xm,
456 ys = m_element_index.ys,
457 ym = m_element_index.ym;
458
459 ParallelSection loop(m_grid->com);
460 try {
461 for (int j=ys; j<ys+ym; j++) {
462 for (int i=xs; i<xs+xm; i++) {
463 // Initialize the map from global to local degrees of freedom for this element.
464 m_element.reset(i, j);
465
466 // Obtain the value of the solution at the nodes adjacent to the element.
467 // Compute the solution values and symmetric gradient at the quadrature points.
468 m_element.nodal_values(*du_local, du_e);
469 if (dirichletBC) {
470 dirichletBC.enforce_homogeneous(m_element, du_e);
471 }
472 quadrature_point_values(m_quadrature, du_e, du_q);
473
474 m_element.nodal_values(u, u_e);
475 if (dirichletBC) {
476 dirichletBC.enforce(m_element, u_e);
477 }
478 quadrature_point_values(m_quadrature, u_e, u_q);
479
480 // Zero out the element-local residual in prep for updating it.
481 for (unsigned int k=0; k<Nk; k++) {
482 dzeta_e[k] = 0;
483 }
484
485 int mask[Nq_max];
486 {
487 Coefficients coeffs[Nk];
488 double thickness[Nq_max];
489 double tauc[Nq_max];
490 double hardness[Nq_max];
491
492 m_element.nodal_values(m_coefficients, coeffs);
493
494 quad_point_values(m_quadrature, coeffs,
495 mask, thickness, tauc, hardness);
496 }
497
498 for (unsigned int q=0; q<Nq; q++) {
499 Vector2 du_qq = du_q[q];
500 Vector2 u_qq = u_q[q];
501
502 // Determine "dbeta/dtauc" at the quadrature point
503 double dbeta_dtauc = 0;
504 if (mask::grounded_ice(mask[q])) {
505 dbeta_dtauc = m_basal_sliding_law->drag(1., u_qq.u, u_qq.v);
506 }
507
508 for (unsigned int k=0; k<Nk; k++) {
509 dzeta_e[k] += W[q]*dbeta_dtauc*(du_qq.u*u_qq.u+du_qq.v*u_qq.v)*test[q][k].val;
510 }
511 } // q
512
513 m_element.add_contribution(dzeta_e, dzeta_a);
514 } // j
515 } // i
516 } catch (...) {
517 loop.failed();
518 }
519 loop.check();
520
521 for (Points p(*m_grid); p; p.next()) {
522 const int i = p.i(), j = p.j();
523
524 double dtauc_dzeta;
525 m_tauc_param.toDesignVariable((*m_zeta)(i, j), NULL, &dtauc_dzeta);
526 dzeta_a[j][i] *= dtauc_dzeta;
527 }
528
529 if (m_fixed_tauc_locations) {
530 fem::DirichletData_Scalar fixedTauc(m_fixed_tauc_locations, NULL);
531 fixedTauc.fix_residual_homogeneous(dzeta_a);
532 }
533 }
534
535 /*!\brief Applies the linearization of the forward map (i.e. the reduced gradient \f$DF\f$ described in
536 the class-level documentation.) */
537 /*! As described previously,
538 \f[
539 Df = J_{\rm State}^{-1} J_{\rm Design}.
540 \f]
541 Applying the linearization then involves the solution of a linear equation.
542 The matrices \f$J_{\rm State}\f$ and \f$J_{\rm Design}\f$ both depend on the value of the
543 design variable \f$\zeta\f$ and the value of the corresponding state variable \f$u=F(\zeta)\f$.
544 These are established by first calling linearize_at.
545 @param[in] dzeta Perturbation of the design variable
546 @param[out] du Computed corresponding perturbation of the state variable; ghosts (if present) are updated.
547 */
548 void IP_SSATaucForwardProblem::apply_linearization(IceModelVec2S &dzeta, IceModelVec2V &du) {
549
550 PetscErrorCode ierr;
551
552 if (m_rebuild_J_state) {
553 this->assemble_jacobian_state(m_velocity, m_J_state);
554 m_rebuild_J_state = false;
555 }
556
557 this->apply_jacobian_design(m_velocity, dzeta, m_du_global);
558 m_du_global.scale(-1);
559
560 // call PETSc to solve linear system by iterative method.
561 ierr = KSPSetOperators(m_ksp, m_J_state, m_J_state);
562 PISM_CHK(ierr, "KSPSetOperators");
563
564 ierr = KSPSolve(m_ksp, m_du_global.vec(), m_du_global.vec());
565 PISM_CHK(ierr, "KSPSolve"); // SOLVE
566
567 KSPConvergedReason reason;
568 ierr = KSPGetConvergedReason(m_ksp, &reason);
569 PISM_CHK(ierr, "KSPGetConvergedReason");
570 if (reason < 0) {
571 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
572 "IP_SSATaucForwardProblem::apply_linearization solve"
573 " failed to converge (KSP reason %s)",
574 KSPConvergedReasons[reason]);
575 } else {
576 m_log->message(4,
577 "IP_SSATaucForwardProblem::apply_linearization converged"
578 " (KSP reason %s)\n",
579 KSPConvergedReasons[reason]);
580 }
581
582 du.copy_from(m_du_global);
583 }
584
585 /*! \brief Applies the transpose of the linearization of the forward map
586 (i.e. the transpose of the reduced gradient \f$DF\f$ described in the class-level documentation.) */
587 /*! As described previously,
588 \f[
589 Df = J_{\rm State}^{-1} J_{\rm Design}.
590 \f]
591 so
592 \f[
593 Df^t = J_{\rm Design}^t \; (J_{\rm State}^t)^{-1} .
594 \f]
595 Applying the transpose of the linearization then involves the solution of a linear equation.
596 The matrices \f$J_{\rm State}\f$ and \f$J_{\rm Design}\f$ both depend on the value of the
597 design variable \f$\zeta\f$ and the value of the corresponding state variable \f$u=F(\zeta)\f$.
598 These are established by first calling linearize_at.
599 @param[in] du Perturbation of the state variable
600 @param[out] dzeta Computed corresponding perturbation of the design variable; ghosts (if present) are updated.
601 */
602 void IP_SSATaucForwardProblem::apply_linearization_transpose(IceModelVec2V &du,
603 IceModelVec2S &dzeta) {
604
605 PetscErrorCode ierr;
606
607 if (m_rebuild_J_state) {
608 this->assemble_jacobian_state(m_velocity, m_J_state);
609 m_rebuild_J_state = false;
610 }
611
612 // Aliases to help with notation consistency below.
613 const IceModelVec2Int *dirichletLocations = m_bc_mask;
614 const IceModelVec2V *dirichletValues = m_bc_values;
615 double dirichletWeight = m_dirichletScale;
616
617 m_du_global.copy_from(du);
618 Vector2 **du_a = m_du_global.get_array();
619 fem::DirichletData_Vector dirichletBC(dirichletLocations, dirichletValues, dirichletWeight);
620
621 if (dirichletBC) {
622 dirichletBC.fix_residual_homogeneous(du_a);
623 }
624
625 m_du_global.end_access();
626
627 // call PETSc to solve linear system by iterative method.
628 ierr = KSPSetOperators(m_ksp, m_J_state, m_J_state);
629 PISM_CHK(ierr, "KSPSetOperators");
630
631 ierr = KSPSolve(m_ksp, m_du_global.vec(), m_du_global.vec());
632 PISM_CHK(ierr, "KSPSolve"); // SOLVE
633
634 KSPConvergedReason reason;
635 ierr = KSPGetConvergedReason(m_ksp, &reason);
636 PISM_CHK(ierr, "KSPGetConvergedReason");
637
638 if (reason < 0) {
639 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
640 "IP_SSATaucForwardProblem::apply_linearization solve"
641 " failed to converge (KSP reason %s)",
642 KSPConvergedReasons[reason]);
643 } else {
644 m_log->message(4,
645 "IP_SSATaucForwardProblem::apply_linearization converged"
646 " (KSP reason %s)\n",
647 KSPConvergedReasons[reason]);
648 }
649
650 this->apply_jacobian_design_transpose(m_velocity, m_du_global, dzeta);
651 dzeta.scale(-1);
652
653 if (dzeta.stencil_width() > 0) {
654 dzeta.update_ghosts();
655 }
656 }
657
658 } // end of namespace inverse
659 } // end of namespace pism