tIP_SSATaucTaoTikhonovProblemLCL.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tIP_SSATaucTaoTikhonovProblemLCL.cc (11456B)
---
1 // Copyright (C) 2012, 2014, 2015, 2016, 2017 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_SSATaucTaoTikhonovProblemLCL.hh"
20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/ConfigInterface.hh"
22
23 namespace pism {
24 namespace inverse {
25
26 typedef IceModelVec2S DesignVec;
27 typedef IceModelVec2V StateVec;
28
29 // typedef TikhonovProblemListener<InverseProblem> Listener;
30 // typedef typename Listener::Ptr ListenerPtr;
31
32 IP_SSATaucTaoTikhonovProblemLCL::IP_SSATaucTaoTikhonovProblemLCL(IP_SSATaucForwardProblem &ssaforward,
33 IP_SSATaucTaoTikhonovProblemLCL::DesignVec &d0,
34 IP_SSATaucTaoTikhonovProblemLCL::StateVec &u_obs,
35 double eta,
36 IPFunctional<DesignVec> &designFunctional,
37 IPFunctional<StateVec> &stateFunctional)
38 : m_ssaforward(ssaforward), m_d0(d0), m_u_obs(u_obs), m_eta(eta),
39 m_designFunctional(designFunctional), m_stateFunctional(stateFunctional) {
40
41 PetscErrorCode ierr;
42 IceGrid::ConstPtr grid = m_d0.grid();
43
44 double stressScale = grid->ctx()->config()->get_number("inverse.design.param_tauc_scale");
45 m_constraintsScale = grid->Lx()*grid->Ly()*4*stressScale;
46
47 m_velocityScale = grid->ctx()->config()->get_number("inverse.ssa.velocity_scale", "m second-1");
48
49
50 int design_stencil_width = m_d0.stencil_width();
51 int state_stencil_width = m_u_obs.stencil_width();
52 m_d.reset(new DesignVec(grid, "design variable", WITH_GHOSTS, design_stencil_width));
53
54 m_d_Jdesign.create(grid, "Jdesign design variable", WITH_GHOSTS, design_stencil_width);
55 m_dGlobal.create(grid, "design variable (global)", WITHOUT_GHOSTS, design_stencil_width);
56 m_dGlobal.copy_from(m_d0);
57
58 m_uGlobal.reset(new StateVec(grid, "state variable (global)",
59 WITHOUT_GHOSTS, state_stencil_width));
60
61 m_u.create(grid, "state variable", WITH_GHOSTS, state_stencil_width);
62 m_du.create(grid, "du", WITH_GHOSTS, state_stencil_width);
63 m_u_Jdesign.create(grid, "Jdesign state variable", WITH_GHOSTS, state_stencil_width);
64
65 m_u_diff.reset(new StateVec(grid, "state residual", WITH_GHOSTS, state_stencil_width));
66
67 m_d_diff.reset(new DesignVec(grid, "design residual", WITH_GHOSTS, design_stencil_width));
68
69 m_dzeta.create(grid,"dzeta",WITH_GHOSTS,design_stencil_width);
70
71 m_grad_state.reset(new StateVec(grid, "state gradient", WITHOUT_GHOSTS, state_stencil_width));
72
73 m_grad_design.reset(new DesignVec(grid, "design gradient", WITHOUT_GHOSTS, design_stencil_width));
74
75 m_constraints.reset(new StateVec(grid,"PDE constraints",WITHOUT_GHOSTS,design_stencil_width));
76
77 DM da;
78 m_ssaforward.get_da(&da);
79
80 ierr = DMSetMatType(da, MATBAIJ);
81 PISM_CHK(ierr, "DMSetMatType");
82
83 ierr = DMCreateMatrix(da, m_Jstate.rawptr());
84 PISM_CHK(ierr, "DMCreateMatrix");
85
86 int nLocalNodes = grid->xm()*grid->ym();
87 int nGlobalNodes = grid->Mx()*grid->My();
88 ierr = MatCreateShell(grid->com, 2*nLocalNodes, nLocalNodes, 2*nGlobalNodes, nGlobalNodes,
89 this, m_Jdesign.rawptr());
90 PISM_CHK(ierr, "MatCreateShell");
91
92 ierr = MatShellSetOperation(m_Jdesign, MATOP_MULT,
93 (void(*)(void))jacobian_design_callback);
94 PISM_CHK(ierr, "MatShellSetOperation");
95
96 ierr = MatShellSetOperation(m_Jdesign, MATOP_MULT_TRANSPOSE,
97 (void(*)(void))jacobian_design_transpose_callback);
98 PISM_CHK(ierr, "MatShellSetOperation");
99
100 m_x.reset(new IPTwoBlockVec(m_dGlobal.vec(),m_uGlobal->vec()));
101 }
102
103 IP_SSATaucTaoTikhonovProblemLCL::~IP_SSATaucTaoTikhonovProblemLCL()
104 {
105 // empty
106 }
107
108 void IP_SSATaucTaoTikhonovProblemLCL::setInitialGuess(DesignVec &d0) {
109 m_dGlobal.copy_from(d0);
110 }
111
112 IP_SSATaucTaoTikhonovProblemLCL::StateVec::Ptr IP_SSATaucTaoTikhonovProblemLCL::stateSolution() {
113
114 m_x->scatterToB(m_uGlobal->vec());
115 m_uGlobal->scale(m_velocityScale);
116
117 return m_uGlobal;
118 }
119
120 IP_SSATaucTaoTikhonovProblemLCL::DesignVec::Ptr IP_SSATaucTaoTikhonovProblemLCL::designSolution() {
121 m_x->scatterToA(m_d->vec()); //CHKERRQ(ierr);
122 return m_d;
123 }
124
125 void IP_SSATaucTaoTikhonovProblemLCL::connect(Tao tao) {
126 PetscErrorCode ierr;
127 ierr = TaoSetStateDesignIS(tao,
128 m_x->blockBIndexSet() /*state*/,
129 m_x->blockAIndexSet() /*design*/);
130 PISM_CHK(ierr, "TaoSetStateDesignIS");
131
132 taoutil::TaoObjGradCallback<IP_SSATaucTaoTikhonovProblemLCL,
133 &IP_SSATaucTaoTikhonovProblemLCL::evaluateObjectiveAndGradient>::connect(tao, *this);
134
135 taoutil::TaoLCLCallbacks<IP_SSATaucTaoTikhonovProblemLCL>::connect(tao, *this,
136 m_constraints->vec(),
137 m_Jstate, m_Jdesign);
138
139 taoutil::TaoMonitorCallback<IP_SSATaucTaoTikhonovProblemLCL>::connect(tao,*this);
140 }
141
142 void IP_SSATaucTaoTikhonovProblemLCL::monitorTao(Tao tao) {
143 PetscErrorCode ierr;
144
145 // Has to be a PetscInt because of the TaoGetSolutionStatus call.
146 PetscInt its;
147 ierr = TaoGetSolutionStatus(tao, &its, NULL, NULL, NULL, NULL, NULL);
148 PISM_CHK(ierr, "TaoGetSolutionStatus");
149
150 int nListeners = m_listeners.size();
151 for (int k = 0; k < nListeners; k++) {
152 m_listeners[k]->iteration(*this, m_eta,
153 its, m_val_design, m_val_state,
154 m_d, m_d_diff, m_grad_design,
155 m_ssaforward.solution(),
156 m_u_diff,
157 m_grad_state,
158 m_constraints);
159 }
160 }
161
162 void IP_SSATaucTaoTikhonovProblemLCL::evaluateObjectiveAndGradient(Tao /*tao*/, Vec x,
163 double *value, Vec gradient) {
164
165 m_x->scatter(x,m_dGlobal.vec(),m_uGlobal->vec());
166 m_uGlobal->scale(m_velocityScale);
167
168 // Variable 'm_dGlobal' has no ghosts. We need ghosts for computation with the design variable.
169 m_d->copy_from(m_dGlobal);
170
171 m_d_diff->copy_from(*m_d);
172 m_d_diff->add(-1,m_d0);
173 m_designFunctional.gradientAt(*m_d_diff, *m_grad_design);
174 m_grad_design->scale(1/m_eta);
175
176 m_u_diff->copy_from(*m_uGlobal);
177 m_u_diff->add(-1, m_u_obs);
178 m_stateFunctional.gradientAt(*m_u_diff, *m_grad_state);
179 m_grad_state->scale(m_velocityScale);
180
181 m_x->gather(m_grad_design->vec(), m_grad_state->vec(), gradient);
182
183 m_designFunctional.valueAt(*m_d_diff, &m_val_design);
184 m_stateFunctional.valueAt(*m_u_diff, &m_val_state);
185
186 *value = m_val_design / m_eta + m_val_state;
187 }
188
189 TerminationReason::Ptr IP_SSATaucTaoTikhonovProblemLCL::formInitialGuess(Vec *x) {
190 m_d->copy_from(m_dGlobal);
191 TerminationReason::Ptr reason = m_ssaforward.linearize_at(*m_d);
192 if (reason->failed()) {
193 return reason;
194 }
195
196 m_uGlobal->copy_from(*m_ssaforward.solution());
197 m_uGlobal->scale(1.0 / m_velocityScale);
198
199 m_x->gather(m_dGlobal.vec(), m_uGlobal->vec());
200
201 // This is probably irrelevant.
202 m_uGlobal->scale(m_velocityScale);
203
204 *x = *m_x;
205 return GenericTerminationReason::success();
206 }
207
208 void IP_SSATaucTaoTikhonovProblemLCL::evaluateConstraints(Tao, Vec x, Vec r) {
209 PetscErrorCode ierr;
210
211 m_x->scatter(x,m_dGlobal.vec(),m_uGlobal->vec());
212 m_uGlobal->scale(m_velocityScale);
213
214 m_d->copy_from(m_dGlobal);
215 m_u.copy_from(*m_uGlobal);
216
217 m_ssaforward.set_design(*m_d);
218
219 m_ssaforward.assemble_residual(m_u, r);
220
221 ierr = VecScale(r,1./m_constraintsScale);
222 PISM_CHK(ierr, "VecScale");
223 }
224
225 void IP_SSATaucTaoTikhonovProblemLCL::evaluateConstraintsJacobianState(Tao, Vec x,
226 Mat Jstate,
227 Mat /*Jpc*/,
228 Mat /*Jinv*/,
229 MatStructure *s) {
230 PetscErrorCode ierr;
231
232 m_x->scatter(x, m_dGlobal.vec(), m_uGlobal->vec());
233 m_uGlobal->scale(m_velocityScale);
234
235 m_d->copy_from(m_dGlobal);
236 m_u.copy_from(*m_uGlobal);
237
238 m_ssaforward.set_design(*m_d);
239
240 m_ssaforward.assemble_jacobian_state(m_u, Jstate);
241
242 *s = SAME_NONZERO_PATTERN;
243
244 ierr = MatScale(Jstate, m_velocityScale / m_constraintsScale);
245 PISM_CHK(ierr, "MatScale");
246 }
247
248 void IP_SSATaucTaoTikhonovProblemLCL::evaluateConstraintsJacobianDesign(Tao, Vec x, Mat /*Jdesign*/) {
249 // I'm not sure if the following are necessary (i.e. will the copies that happen
250 // in evaluateObjectiveAndGradient be sufficient) but we'll do them here
251 // just in case.
252 m_x->scatter(x,m_dGlobal.vec(),m_uGlobal->vec());
253 m_uGlobal->scale(m_velocityScale);
254 m_d_Jdesign.copy_from(m_dGlobal);
255 m_u_Jdesign.copy_from(*m_uGlobal);
256 }
257
258 void IP_SSATaucTaoTikhonovProblemLCL::applyConstraintsJacobianDesign(Vec x, Vec y) {
259 m_dzeta.copy_from_vec(x);
260
261 m_ssaforward.set_design(m_d_Jdesign);
262
263 m_ssaforward.apply_jacobian_design(m_u_Jdesign, m_dzeta, y);
264
265 PetscErrorCode ierr = VecScale(y,1./m_constraintsScale);
266 PISM_CHK(ierr, "VecScale");
267 }
268
269 void IP_SSATaucTaoTikhonovProblemLCL::applyConstraintsJacobianDesignTranspose(Vec x, Vec y) {
270 m_du.copy_from_vec(x);
271
272 m_ssaforward.set_design(m_d_Jdesign);
273
274 m_ssaforward.apply_jacobian_design_transpose(m_u_Jdesign, m_du, y);
275
276 PetscErrorCode ierr = VecScale(y, 1.0 / m_constraintsScale);
277 PISM_CHK(ierr, "VecScale");
278 }
279
280 PetscErrorCode IP_SSATaucTaoTikhonovProblemLCL::jacobian_design_callback(Mat A, Vec x, Vec y) {
281 try {
282 IP_SSATaucTaoTikhonovProblemLCL *ctx;
283 PetscErrorCode ierr = MatShellGetContext(A,&ctx);
284 PISM_CHK(ierr, "MatShellGetContext");
285
286 ctx->applyConstraintsJacobianDesign(x,y);
287 } catch (...) {
288 MPI_Comm com = MPI_COMM_SELF;
289 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
290 handle_fatal_errors(com);
291 SETERRQ(com, 1, "A PISM callback failed");
292 }
293 return 0;
294 }
295
296 PetscErrorCode IP_SSATaucTaoTikhonovProblemLCL::jacobian_design_transpose_callback(Mat A, Vec x,
297 Vec y) {
298 try {
299 IP_SSATaucTaoTikhonovProblemLCL *ctx;
300 PetscErrorCode ierr = MatShellGetContext(A,&ctx);
301 PISM_CHK(ierr, "MatShellGetContext");
302
303 ctx->applyConstraintsJacobianDesignTranspose(x,y);
304 } catch (...) {
305 MPI_Comm com = MPI_COMM_SELF;
306 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
307 handle_fatal_errors(com);
308 SETERRQ(com, 1, "A PISM callback failed");
309 }
310 return 0;
311 }
312
313 } // end of namespace inverse
314 } // end of namespace pism