PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
IP_SSATaucTaoTikhonovProblemLCL.cc
Go to the documentation of this file.
1 // Copyright (C) 2012, 2014, 2015, 2016, 2017, 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/util/IceGrid.hh"
21 #include "pism/util/ConfigInterface.hh"
22 #include "pism/util/Context.hh"
23 
24 namespace pism {
25 namespace inverse {
26 
29 
30 // typedef TikhonovProblemListener<InverseProblem> Listener;
31 // typedef typename Listener::Ptr ListenerPtr;
32 
36  double eta,
37  IPFunctional<DesignVec> &designFunctional,
38  IPFunctional<StateVec> &stateFunctional)
39 : m_ssaforward(ssaforward),
40  m_dGlobal(d0.grid(), "design variable (global)", WITHOUT_GHOSTS, d0.stencil_width()),
41  m_d0(d0),
42  m_dzeta(d0.grid(),"dzeta",WITH_GHOSTS, d0.stencil_width()),
43  m_u(d0.grid(), "state variable", WITH_GHOSTS, u_obs.stencil_width()),
44  m_du(d0.grid(), "du", WITH_GHOSTS, u_obs.stencil_width()),
45  m_u_obs(u_obs),
46  m_eta(eta),
47  m_d_Jdesign(d0.grid(), "Jdesign design variable", WITH_GHOSTS, d0.stencil_width()),
48  m_u_Jdesign(d0.grid(), "Jdesign state variable", WITH_GHOSTS, u_obs.stencil_width()),
49  m_designFunctional(designFunctional),
50  m_stateFunctional(stateFunctional)
51 {
52 
53  PetscErrorCode ierr;
54  IceGrid::ConstPtr grid = m_d0.grid();
55 
56  double stressScale = grid->ctx()->config()->get_number("inverse.design.param_tauc_scale");
57  m_constraintsScale = grid->Lx()*grid->Ly()*4*stressScale;
58 
59  m_velocityScale = grid->ctx()->config()->get_number("inverse.ssa.velocity_scale", "m second-1");
60 
61 
62  int design_stencil_width = m_d0.stencil_width();
63  int state_stencil_width = m_u_obs.stencil_width();
64  m_d.reset(new DesignVec(grid, "design variable", WITH_GHOSTS, design_stencil_width));
65 
67 
68  m_uGlobal.reset(new StateVec(grid, "state variable (global)",
69  WITHOUT_GHOSTS, state_stencil_width));
70 
71  m_u_diff.reset(new StateVec(grid, "state residual", WITH_GHOSTS, state_stencil_width));
72 
73  m_d_diff.reset(new DesignVec(grid, "design residual", WITH_GHOSTS, design_stencil_width));
74 
75  m_grad_state.reset(new StateVec(grid, "state gradient", WITHOUT_GHOSTS, state_stencil_width));
76 
77  m_grad_design.reset(new DesignVec(grid, "design gradient", WITHOUT_GHOSTS, design_stencil_width));
78 
79  m_constraints.reset(new StateVec(grid,"PDE constraints",WITHOUT_GHOSTS,design_stencil_width));
80 
82 
83  ierr = DMSetMatType(da, MATBAIJ);
84  PISM_CHK(ierr, "DMSetMatType");
85 
86  ierr = DMCreateMatrix(da, m_Jstate.rawptr());
87  PISM_CHK(ierr, "DMCreateMatrix");
88 
89  int nLocalNodes = grid->xm()*grid->ym();
90  int nGlobalNodes = grid->Mx()*grid->My();
91  ierr = MatCreateShell(grid->com, 2*nLocalNodes, nLocalNodes, 2*nGlobalNodes, nGlobalNodes,
92  this, m_Jdesign.rawptr());
93  PISM_CHK(ierr, "MatCreateShell");
94 
95  ierr = MatShellSetOperation(m_Jdesign, MATOP_MULT,
96  (void(*)(void))jacobian_design_callback);
97  PISM_CHK(ierr, "MatShellSetOperation");
98 
99  ierr = MatShellSetOperation(m_Jdesign, MATOP_MULT_TRANSPOSE,
100  (void(*)(void))jacobian_design_transpose_callback);
101  PISM_CHK(ierr, "MatShellSetOperation");
102 
103  m_x.reset(new IPTwoBlockVec(m_dGlobal.vec(),m_uGlobal->vec()));
104 }
105 
107  m_dGlobal.copy_from(d0);
108 }
109 
111 
112  m_x->scatterToB(m_uGlobal->vec());
113  m_uGlobal->scale(m_velocityScale);
114 
115  return m_uGlobal;
116 }
117 
119  m_x->scatterToA(m_d->vec()); //CHKERRQ(ierr);
120  return m_d;
121 }
122 
124  PetscErrorCode ierr;
125  ierr = TaoSetStateDesignIS(tao,
126  m_x->blockBIndexSet() /*state*/,
127  m_x->blockAIndexSet() /*design*/);
128  PISM_CHK(ierr, "TaoSetStateDesignIS");
129 
132 
134  m_constraints->vec(),
136 
138 }
139 
141  PetscErrorCode ierr;
142 
143  // Has to be a PetscInt because of the TaoGetSolutionStatus call.
144  PetscInt its;
145  ierr = TaoGetSolutionStatus(tao, &its, NULL, NULL, NULL, NULL, NULL);
146  PISM_CHK(ierr, "TaoGetSolutionStatus");
147 
148  int nListeners = m_listeners.size();
149  for (int k = 0; k < nListeners; k++) {
150  m_listeners[k]->iteration(*this, m_eta,
154  m_u_diff,
155  m_grad_state,
156  m_constraints);
157  }
158 }
159 
161  double *value, Vec gradient) {
162 
163  m_x->scatter(x,m_dGlobal.vec(),m_uGlobal->vec());
164  m_uGlobal->scale(m_velocityScale);
165 
166  // Variable 'm_dGlobal' has no ghosts. We need ghosts for computation with the design variable.
167  m_d->copy_from(m_dGlobal);
168 
169  m_d_diff->copy_from(*m_d);
170  m_d_diff->add(-1,m_d0);
172  m_grad_design->scale(1/m_eta);
173 
174  m_u_diff->copy_from(*m_uGlobal);
175  m_u_diff->add(-1, m_u_obs);
178 
179  m_x->gather(m_grad_design->vec(), m_grad_state->vec(), gradient);
180 
183 
184  *value = m_val_design / m_eta + m_val_state;
185 }
186 
188  m_d->copy_from(m_dGlobal);
190  if (reason->failed()) {
191  return reason;
192  }
193 
194  m_uGlobal->copy_from(*m_ssaforward.solution());
195  m_uGlobal->scale(1.0 / m_velocityScale);
196 
197  m_x->gather(m_dGlobal.vec(), m_uGlobal->vec());
198 
199  // This is probably irrelevant.
200  m_uGlobal->scale(m_velocityScale);
201 
202  *x = *m_x;
204 }
205 
207  PetscErrorCode ierr;
208 
209  m_x->scatter(x,m_dGlobal.vec(),m_uGlobal->vec());
210  m_uGlobal->scale(m_velocityScale);
211 
212  m_d->copy_from(m_dGlobal);
214 
216 
218 
219  ierr = VecScale(r,1./m_constraintsScale);
220  PISM_CHK(ierr, "VecScale");
221 }
222 
224  Mat Jstate,
225  Mat /*Jpc*/,
226  Mat /*Jinv*/,
227  MatStructure *s) {
228  PetscErrorCode ierr;
229 
230  m_x->scatter(x, m_dGlobal.vec(), m_uGlobal->vec());
231  m_uGlobal->scale(m_velocityScale);
232 
233  m_d->copy_from(m_dGlobal);
235 
237 
239 
240  *s = SAME_NONZERO_PATTERN;
241 
242  ierr = MatScale(Jstate, m_velocityScale / m_constraintsScale);
243  PISM_CHK(ierr, "MatScale");
244 }
245 
247  // I'm not sure if the following are necessary (i.e. will the copies that happen
248  // in evaluateObjectiveAndGradient be sufficient) but we'll do them here
249  // just in case.
250  m_x->scatter(x,m_dGlobal.vec(),m_uGlobal->vec());
251  m_uGlobal->scale(m_velocityScale);
254 }
255 
257 
258  PetscErrorCode ierr;
259  {
260  ierr = DMGlobalToLocalBegin(*m_dzeta.dm(), x, INSERT_VALUES, m_dzeta.vec());
261  PISM_CHK(ierr, "DMGlobalToLocalBegin");
262 
263  ierr = DMGlobalToLocalEnd(*m_dzeta.dm(), x, INSERT_VALUES, m_dzeta.vec());
264  PISM_CHK(ierr, "DMGlobalToLocalEnd");
265  }
266 
268 
270 
271  ierr = VecScale(y, 1.0 / m_constraintsScale); PISM_CHK(ierr, "VecScale");
272 }
273 
275 
276  PetscErrorCode ierr;
277  {
278  ierr = DMGlobalToLocalBegin(*m_du.dm(), x, INSERT_VALUES, m_du.vec());
279  PISM_CHK(ierr, "DMGlobalToLocalBegin");
280 
281  ierr = DMGlobalToLocalEnd(*m_du.dm(), x, INSERT_VALUES, m_du.vec());
282  PISM_CHK(ierr, "DMGlobalToLocalEnd");
283  }
284 
286 
288 
289  ierr = VecScale(y, 1.0 / m_constraintsScale); PISM_CHK(ierr, "VecScale");
290 }
291 
293  try {
295  PetscErrorCode ierr = MatShellGetContext(A,&ctx);
296  PISM_CHK(ierr, "MatShellGetContext");
297 
299  } catch (...) {
300  MPI_Comm com = MPI_COMM_SELF;
301  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
302  handle_fatal_errors(com);
303  SETERRQ(com, 1, "A PISM callback failed");
304  }
305  return 0;
306 }
307 
309  Vec y) {
310  try {
312  PetscErrorCode ierr = MatShellGetContext(A,&ctx);
313  PISM_CHK(ierr, "MatShellGetContext");
314 
316  } catch (...) {
317  MPI_Comm com = MPI_COMM_SELF;
318  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
319  handle_fatal_errors(com);
320  SETERRQ(com, 1, "A PISM callback failed");
321  }
322  return 0;
323 }
324 
325 } // end of namespace inverse
326 } // end of namespace pism
static TerminationReason::Ptr success()
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
void copy_from(const IceModelVec2S &source)
std::shared_ptr< IceModelVec2S > Ptr
Definition: iceModelVec.hh:341
std::shared_ptr< IceModelVec2V > Ptr
void copy_from(const IceModelVec2< T > &source)
Definition: IceModelVec2.hh:97
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
Definition: iceModelVec.cc:334
petsc::Vec & vec() const
Definition: iceModelVec.cc:342
std::shared_ptr< petsc::DM > dm() const
Definition: iceModelVec.cc:356
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
std::shared_ptr< TerminationReason > Ptr
virtual void valueAt(IMVecType &x, double *OUTPUT)=0
Computes the value of the functional at the vector x.
virtual void gradientAt(IMVecType &x, IMVecType &gradient)=0
Computes the gradient of the functional at the vector x.
Abstract base class for functions from ice model vectors to .
Definition: IPFunctional.hh:40
virtual void set_design(IceModelVec2S &zeta)
Sets the current value of of the design paramter .
virtual void assemble_residual(IceModelVec2V &u, IceModelVec2V &R)
Computes the residual function as defined in the class-level documentation.
virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, IceModelVec2V &du)
Applies the design Jacobian matrix to a perturbation of the design variable.
petsc::DM & get_da() const
Exposes the DMDA of the underlying grid for the benefit of TAO.
virtual IceModelVec2V::Ptr solution()
Returns the last solution of the SSA as computed by linearize_at.
virtual void assemble_jacobian_state(IceModelVec2V &u, Mat J)
Assembles the state Jacobian matrix.
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.
virtual TerminationReason::Ptr linearize_at(IceModelVec2S &zeta)
Sets the current value of the design variable and solves the SSA to find the associated .
Implements the forward problem of the map taking to the corresponding solution of the SSA.
virtual void evaluateObjectiveAndGradient(Tao tao, Vec x, double *value, Vec gradient)
static PetscErrorCode jacobian_design_transpose_callback(Mat A, Vec x, Vec y)
virtual void evaluateConstraintsJacobianState(Tao tao, Vec x, Mat Jstate, Mat Jpc, Mat Jinv, MatStructure *s)
IP_SSATaucTaoTikhonovProblemLCL(IP_SSATaucForwardProblem &ssaforward, DesignVec &d0, StateVec &u_obs, double eta, IPFunctional< DesignVec > &designFunctional, IPFunctional< StateVec > &stateFunctional)
virtual void evaluateConstraintsJacobianDesign(Tao tao, Vec x, Mat Jdesign)
static PetscErrorCode jacobian_design_callback(Mat A, Vec x, Vec y)
Defines a Tikhonov minimization problem of determining from SSA velocities to be solved with a TaoBa...
static void connect(Tao tao, Problem &p, Vec c, Mat Jc, Mat Jd, Mat Jcpc=NULL, Mat Jcinv=NULL)
Definition: TaoUtil.hh:459
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:232
Adaptor to connect a TAO objective and gradient function callback to a C++ object method.
Definition: TaoUtil.hh:407
#define PISM_CHK(errcode, name)
static const double k
Definition: exactTestP.cc:45
void handle_fatal_errors(MPI_Comm com)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49