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