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