PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
IPTaoTikhonovProblem.hh
Go to the documentation of this file.
1 // Copyright (C) 2012,2013,2014,2015,2016,2017,2020,2022,2023 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 #ifndef PISM_IPTAOTIKHONOVPROBLEM_HH
20 #define PISM_IPTAOTIKHONOVPROBLEM_HH
21 
22 #include <memory>
23 
24 #include "pism/inverse/TaoUtil.hh"
25 #include "pism/inverse/functional/IPFunctional.hh"
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/util/Context.hh"
28 #include "pism/util/Grid.hh"
29 #include "pism/util/Logger.hh"
30 #include "pism/util/array/Array.hh" //
31 #include "pism/util/petscwrappers/DM.hh"
32 #include "pism/util/petscwrappers/Vec.hh"
33 
34 namespace pism {
35 
36 class Grid;
37 
38 namespace inverse {
39 
40 template <class ForwardProblem>
41 class IPTaoTikhonovProblem;
42 
43 //! Iteration callback class for IPTaoTikhonovProblem
44 /** A class for objects receiving iteration callbacks from a
45  * IPTaoTikhonovProblem. These callbacks can be used to monitor the
46  * solution, plot iterations, print diagnostic messages, etc.
47  * IPTaoTikhonovProblemListeners are ususally used via a reference
48  * counted pointer IPTaoTikhonovProblemListener::Ptr to allow for good
49  * memory management when Listeners are created as subclasses of
50  * python classes. It would have been better to nest this inside of
51  * IPTaoTikhonovProblem, but SWIG has a hard time with nested classes,
52  * so it's outer instead.
53  */
54 template <class ForwardProblem>
56 public:
57  typedef std::shared_ptr<IPTaoTikhonovProblemListener> Ptr;
58 
59  typedef std::shared_ptr<typename ForwardProblem::DesignVec> DesignVecPtr;
60  typedef std::shared_ptr<typename ForwardProblem::StateVec> StateVecPtr;
61  typedef std::shared_ptr<typename ForwardProblem::StateVec1> StateVec1Ptr;
62 
64  }
66  }
67 
68  //! The method called after each minimization iteration.
69  virtual void iteration(IPTaoTikhonovProblem<ForwardProblem> &problem, double eta, int iter,
70  double objectiveValue, double designValue, const DesignVecPtr d,
71  const DesignVecPtr diff_d, const DesignVecPtr grad_d, const StateVecPtr u,
72  const StateVecPtr diff_u, const DesignVecPtr grad_u,
73  const DesignVecPtr gradient) = 0;
74 };
75 
76 
77 //! \brief Defines a Tikhonov minimization problem to be solved with a TaoBasicSolver.
78 /*! Suppose \f$F\f$ is a map from a space \f$D\f$ of design variables to a space \f$S\f$ of
79  state variables and we wish to solve a possibly ill-posed problem of the form
80  \f[ F(d) = u \f]
81  where \f$u\f$ is know and \f$d\f$ is unknown. Approximate solutions can be obtained by
82  finding minimizers of an associated Tikhonov functional
83  \f[
84  J(d) = J_{S}(F(d)-u) + \frac{1}{\eta}J_{D}(d-d_0)
85  \f]
86  where \$J_{D}\$ and \$J_{S}\$ are functionals on the spaces \f$D\f$ and \f$S\f$ respectively,
87  \f$\eta\f$ is a penalty parameter, and \f$d_0\f$ is a best a-priori guess for the the solution.
88  The IPTaoTikhonovProblem class encapuslates all of the data required to formulate the minimization
89  problem as a Problem tha can be solved using a TaoBasicSolver. It is templated on the
90  the class ForwardProblem which defines the class of the forward map \f$F\f$ as well as the
91  spaces \f$D\f$ and \f$S\f$. An instance of ForwardProblem, along with
92  specific functionals \f$J_D\f$ and \f$J_S\f$, the parameter \f$\eta\f$, and the data
93  \f$y\f$ and \f$x_0\f$ are provided on constructing a IPTaoTikhonovProblem.
94 
95  For example, if the SSATaucForwardProblem class defines the map taking yield stresses \f$\tau_c\f$
96  to the corresponding surface velocity field solving the SSA, a schematic setup of solving
97  the associated Tikhonov problem goes as follows.
98 
99  \code
100  SSATaucForwardProblem forwardProblem(grid);
101  L2NormFunctional2S designFunctional(grid); //J_X
102  L2NormFunctional2V stateFunctional(grid); //J_Y
103  array::Vector u_obs; // Set this to the surface velocity observations.
104  array::Scalar tauc_0; // Set this to the initial guess for tauc.
105  double eta; // Set this to the desired penalty parameter.
106 
107  typedef InvSSATauc IPTaoTikhonovProblem<SSATaucForwardProblem>;
108  InvSSATauc tikhonovProblem(forwardProblem,tauc_0,u_obs,eta,designFunctional,stateFunctional);
109 
110  TaoBasicSolver<InvSSATauc> solver(com, "cg", tikhonovProblem);
111 
112  TerminationReason::Ptr reason = solver.solve();
113 
114  if (reason->succeeded()) {
115  printf("Success: %s\n",reason->description().c_str());
116  } else {
117  printf("Failure: %s\n",reason->description().c_str());
118  }
119  \endcode
120 
121  The class ForwardProblem that defines the forward problem
122  must have the following characteristics:
123 
124  <ol>
125  <li> Contains typedefs for DesignVec and StateVec that effectively
126  define the function spaces \f$D\f$ and \f$S\f$. E.g.
127 
128  \code
129  typedef array::Scalar DesignVec;
130  typedef array::Vector StateVec;
131  \endcode
132  would be appropriate for a map from basal yeild stress to surface velocities.
133 
134  <li> A method
135  \code
136  TerminationReason::Ptr linearize_at(DesignVec &d);
137  \endcode
138  that instructs the class to compute the value of F and
139  anything needed to compute its linearization at \a d. This is the first method
140  called when working with a new iterate of \a d.
141 
142  <li> A method
143  \code
144  StateVec &solution()
145  \endcode
146  that returns the most recently computed value of \f$F(d)\f$
147  as computed by a call to linearize_at.
148 
149  <li> A method
150  \code
151  void apply_linearization_transpose(StateVec &du, DesignVec &dzeta);
152  \endcode
153  that computes the action of \f$(F')^t\f$,
154  where \f$F'\f$ is the linearization of \f$F\f$ at the current iterate, and the transpose
155  is computed in the standard sense (i.e. thinking of \f$F'\f$ as a matrix with respect
156  to the bases implied by the DesignVec and StateVec spaces). The need for a transpose arises
157  because
158  \f[
159  \frac{d}{dt} J_{S}(F(d+t\delta d)-u) = [DJ_S]_{k}\; F'_{kj} \; \delta d
160  \f]
161  and hence the gradient of the term \f$J_{S}(F(d)-u)\f$ with respect to \f$d\f$ is given
162  by
163  \f[
164  (F')^t (\nabla J_S)^t.
165  \f]
166  </ol>
167 */
168 template <class ForwardProblem>
170 public:
173  typedef typename ForwardProblem::StateVec1 StateVec1;
174 
175  typedef typename ForwardProblem::DesignVecGhosted DesignVecGhosted;
176  typedef std::shared_ptr<typename ForwardProblem::DesignVecGhosted> DesignVecGhostedPtr;
177 
178  typedef std::shared_ptr<typename ForwardProblem::DesignVec> DesignVecPtr;
179  typedef std::shared_ptr<typename ForwardProblem::StateVec> StateVecPtr;
180  typedef std::shared_ptr<typename ForwardProblem::StateVec1> StateVec1Ptr;
181 
182  /*! Constructs a Tikhonov problem:
183 
184  Minimize \f$J(d) = J_S(F(d)-u_obs) + \frac{1}{\eta} J_D(d-d0) \f$
185 
186  that can be solved with a TaoBasicSolver.
187 
188  @param forward Class defining the map F. See class-level documentation for requirements of F.
189  @param d0 Best a-priori guess for the design parameter.
190  @param u_obs State parameter to match (i.e. approximately solve F(d)=u_obs)
191  @param eta Penalty parameter/Lagrange multiplier. Take eta to zero to impose more regularization to an ill posed problem.
192  @param designFunctional The functional \f$J_D\f$
193  @param stateFunctional The functional \f$J_S\f$
194  */
195 
196  IPTaoTikhonovProblem(ForwardProblem &forward, DesignVec &d0, StateVec &u_obs, double eta,
197  IPFunctional<DesignVec> &designFunctional,
198  IPFunctional<StateVec> &stateFunctional);
199 
201 
202 
203  //! Sets the initial guess for minimization iterations. If this isn't set explicitly,
204  // the parameter \f$d0\f$ appearing the in the Tikhonov functional will be used.
205  virtual void setInitialGuess(DesignVec &d) {
206  m_dGlobal.copy_from(d);
207  }
208 
209  //! Callback provided to TAO for objective evaluation.
210  virtual void evaluateObjectiveAndGradient(Tao tao, Vec x, double *value, Vec gradient);
211 
212  //! Add an object to the list of objects to be called after each iteration.
214  m_listeners.push_back(listener);
215  }
216 
217  //! Final value of \f$F(d)\f$, where \f$d\f$ is the solution of the minimization.
219  return m_forward.solution();
220  }
221 
222  //! Value of \f$d\f$, the solution of the minimization problem.
224  return m_d;
225  }
226 
227  //! Callback from TaoBasicSolver, used to wire the connections between a Tao and
228  // the current class.
229  virtual void connect(Tao tao);
230 
231  //! Callback from TAO after each iteration. The call is forwarded to each element of our list of listeners.
232  virtual void monitorTao(Tao tao);
233 
234  //! Callback from TAO to detect convergence. Allows us to implement a custom convergence check.
235  virtual void convergenceTest(Tao tao);
236 
237  //! Callback from TaoBasicSolver to form the starting iterate for the minimization. See also
238  // setInitialGuess.
239  virtual std::shared_ptr<TerminationReason> formInitialGuess(Vec *v) {
240  *v = m_dGlobal.vec();
242  }
243 
244 protected:
245  std::shared_ptr<const Grid> m_grid;
246 
247  ForwardProblem &m_forward;
248 
249  /// Current iterate of design parameter
251  /// Initial iterate of design parameter, stored without ghosts for the benefit of TAO.
253  /// A-priori estimate of design parameter
255  /// Storage for (m_d-m_d0)
256  DesignVecPtr m_d_diff; // ghosted
257 
258  /// State parameter to match via F(d)=u_obs
260  /// Storage for F(d)-u_obs
261  StateVec1Ptr m_u_diff; // ghosted
262 
263  /// Temporary storage used in gradient computation.
265 
266  /// Gradient of \f$J_D\f$ at the current iterate.
268  /// Gradient of \f$J_S\f$ at the current iterate.
270  /** Weighted sum of the design and state gradients corresponding to
271  * the gradient of the Tikhonov functional \f$J\f$.
272  */
274 
275  /// Penalty parameter/Lagrange multiplier.
276  double m_eta;
277 
278  /// Value of \f$J_D\f$ at the current iterate.
279  double m_val_design;
280  /// Value of \f$J_S\f$ at the current iterate.
281  double m_val_state;
282 
283  /// Implementation of \f$J_D\f$.
285  /// Implementation of \f$J_S\f$.
287 
288  /// List of iteration callbacks.
289  std::vector<typename IPTaoTikhonovProblemListener<ForwardProblem>::Ptr> m_listeners;
290 
291  /// Convergence parameter: convergence stops when \f$||J_D||_2 <\f$ m_tikhonov_rtol.
293 
294  /** Convergence parameter: convergence stops when \f$||J_D||_2 \f$
295  * is less than m_tikhonov_rtol times the maximum of the gradient of
296  * \f$J_S\f$ and \f$(1/\eta)J_D\f$. This occurs when the two terms
297  * forming the sum of the gradient of \f$J\f$ point in roughly
298  * opposite directions with the same magnitude.
299  */
301 };
302 
303 template <class ForwardProblem>
305  ForwardProblem &forward, DesignVec &d0, StateVec &u_obs, double eta,
306  IPFunctional<DesignVec> &designFunctional, IPFunctional<StateVec> &stateFunctional)
307  : m_grid(d0.grid()),
308  m_forward(forward),
309  m_dGlobal(d0.grid(), "design variable (global)"),
310  m_d0(d0),
311  m_u_obs(u_obs),
312  m_adjointRHS(d0.grid(), "work vector"),
313  m_eta(eta),
314  m_designFunctional(designFunctional),
315  m_stateFunctional(stateFunctional) {
316 
317  m_tikhonov_atol = m_grid->ctx()->config()->get_number("inverse.tikhonov.atol");
318  m_tikhonov_rtol = m_grid->ctx()->config()->get_number("inverse.tikhonov.rtol");
319 
320  m_d = std::make_shared<DesignVecGhosted>(m_grid, "design variable");
321 
322  m_dGlobal.copy_from(m_d0);
323 
324  m_u_diff = std::make_shared<StateVec1>(m_grid, "state residual");
325 
326  m_d_diff = std::make_shared<DesignVecGhosted>(m_grid, "design residual");
327 
328  m_grad_state = std::make_shared<DesignVec>(m_grid, "state gradient");
329 
330  m_grad_design = std::make_shared<DesignVec>(m_grid, "design gradient");
331 
332  m_grad = std::make_shared<DesignVec>(m_grid, "gradient");
333 }
334 
335 template<class ForwardProblem>
337  // empty
338 }
339 
340 template<class ForwardProblem>
344 
345  ObjGradCallback::connect(tao,*this);
346 
348 
350 
351  double
352  gatol = PETSC_DEFAULT,
353  grtol = PETSC_DEFAULT,
354  gttol = PETSC_DEFAULT;
355 
356 #if PETSC_VERSION_LT(3,7,0)
357  double fatol = 1e-10, frtol = 1e-20;
358  PetscErrorCode ierr = TaoSetTolerances(tao, fatol, frtol, gatol, grtol, gttol);
359  PISM_CHK(ierr, "TaoSetTolerances");
360 #else
361  PetscErrorCode ierr = TaoSetTolerances(tao, gatol, grtol, gttol);
362  PISM_CHK(ierr, "TaoSetTolerances");
363 #endif
364 }
365 
366 template<class ForwardProblem>
368  PetscInt its;
369  TaoGetSolutionStatus(tao, &its, NULL, NULL, NULL, NULL, NULL);
370 
371  int nListeners = m_listeners.size();
372  for (int k=0; k<nListeners; k++) {
373  m_listeners[k]->iteration(*this, m_eta,
374  its, m_val_design, m_val_state,
375  m_d, m_d_diff, m_grad_design,
376  m_forward.solution(), m_u_diff, m_grad_state,
377  m_grad);
378  }
379 }
380 
381 template<class ForwardProblem> void IPTaoTikhonovProblem<ForwardProblem>::convergenceTest(Tao tao) {
382  double designNorm, stateNorm, sumNorm;
383  double dWeight, sWeight;
384  dWeight = 1/m_eta;
385  sWeight = 1;
386 
387  designNorm = m_grad_design->norm(NORM_2)[0];
388  stateNorm = m_grad_state->norm(NORM_2)[0];
389  sumNorm = m_grad->norm(NORM_2)[0];
390 
391  designNorm *= dWeight;
392  stateNorm *= sWeight;
393 
394  if (sumNorm < m_tikhonov_atol) {
395  TaoSetConvergedReason(tao, TAO_CONVERGED_GATOL);
396  } else if (sumNorm < m_tikhonov_rtol*std::max(designNorm,stateNorm)) {
397  TaoSetConvergedReason(tao, TAO_CONVERGED_USER);
398  } else {
399  TaoDefaultConvergenceTest(tao, NULL);
400  }
401 }
402 
403 template<class ForwardProblem>
405  double *value, Vec gradient) {
406  PetscErrorCode ierr;
407  // Variable 'x' has no ghosts. We need ghosts for computation with the design variable.
408  {
409  ierr = DMGlobalToLocalBegin(*m_d->dm(), x, INSERT_VALUES, m_d->vec());
410  PISM_CHK(ierr, "DMGlobalToLocalBegin");
411 
412  ierr = DMGlobalToLocalEnd(*m_d->dm(), x, INSERT_VALUES, m_d->vec());
413  PISM_CHK(ierr, "DMGlobalToLocalEnd");
414  }
415 
416  auto reason = m_forward.linearize_at(*m_d);
417  if (reason->failed()) {
418  Logger::ConstPtr log = m_grid->ctx()->log();
419  log->message(2,
420  "IPTaoTikhonovProblem::evaluateObjectiveAndGradient"
421  " failure in forward solve\n%s\n", reason->description().c_str());
422  ierr = TaoSetConvergedReason(tao, TAO_DIVERGED_USER);
423  PISM_CHK(ierr, "TaoSetConvergedReason");
424  return;
425  }
426 
427  m_d_diff->copy_from(*m_d);
428  m_d_diff->add(-1, m_d0);
429  m_designFunctional.gradientAt(*m_d_diff, *m_grad_design);
430 
431  m_u_diff->copy_from(*m_forward.solution());
432  m_u_diff->add(-1, m_u_obs);
433 
434  // The following computes the reduced gradient.
435  m_stateFunctional.gradientAt(*m_u_diff, m_adjointRHS);
436  m_forward.apply_linearization_transpose(m_adjointRHS, *m_grad_state);
437 
438  m_grad->copy_from(*m_grad_design);
439  m_grad->scale(1.0 / m_eta);
440  m_grad->add(1, *m_grad_state);
441 
442  ierr = VecCopy(m_grad->vec(), gradient);
443  PISM_CHK(ierr, "VecCopy");
444 
445  double valDesign, valState;
446  m_designFunctional.valueAt(*m_d_diff, &valDesign);
447  m_stateFunctional.valueAt(*m_u_diff, &valState);
448 
449  m_val_design = valDesign;
450  m_val_state = valState;
451 
452  *value = valDesign / m_eta + valState;
453 }
454 
455 } // end of namespace inverse
456 } // end of namespace pism
457 
458 #endif // PISM_IPTAOTIKHONOVPROBLEM_HH
static std::shared_ptr< TerminationReason > success()
std::shared_ptr< const Logger > ConstPtr
Definition: Logger.hh:46
Abstract base class for functions from ice model vectors to .
Definition: IPFunctional.hh:41
std::shared_ptr< IPTaoTikhonovProblemListener > Ptr
virtual void iteration(IPTaoTikhonovProblem< ForwardProblem > &problem, double eta, int iter, double objectiveValue, double designValue, const DesignVecPtr d, const DesignVecPtr diff_d, const DesignVecPtr grad_d, const StateVecPtr u, const StateVecPtr diff_u, const DesignVecPtr grad_u, const DesignVecPtr gradient)=0
The method called after each minimization iteration.
std::shared_ptr< typename ForwardProblem::StateVec1 > StateVec1Ptr
std::shared_ptr< typename ForwardProblem::StateVec > StateVecPtr
std::shared_ptr< typename ForwardProblem::DesignVec > DesignVecPtr
Iteration callback class for IPTaoTikhonovProblem.
double m_eta
Penalty parameter/Lagrange multiplier.
std::shared_ptr< typename ForwardProblem::DesignVec > DesignVecPtr
std::vector< typename IPTaoTikhonovProblemListener< ForwardProblem >::Ptr > m_listeners
List of iteration callbacks.
virtual void convergenceTest(Tao tao)
Callback from TAO to detect convergence. Allows us to implement a custom convergence check.
std::shared_ptr< typename ForwardProblem::StateVec1 > StateVec1Ptr
double m_val_state
Value of at the current iterate.
virtual void addListener(typename IPTaoTikhonovProblemListener< ForwardProblem >::Ptr listener)
Add an object to the list of objects to be called after each iteration.
double m_tikhonov_atol
Convergence parameter: convergence stops when m_tikhonov_rtol.
DesignVecPtr m_grad_state
Gradient of at the current iterate.
virtual void setInitialGuess(DesignVec &d)
Sets the initial guess for minimization iterations. If this isn't set explicitly,.
DesignVecPtr m_d_diff
Storage for (m_d-m_d0)
virtual void connect(Tao tao)
Callback from TaoBasicSolver, used to wire the connections between a Tao and.
std::shared_ptr< typename ForwardProblem::StateVec > StateVecPtr
IPFunctional< array::Scalar > & m_designFunctional
Implementation of .
virtual void evaluateObjectiveAndGradient(Tao tao, Vec x, double *value, Vec gradient)
Callback provided to TAO for objective evaluation.
StateVec1Ptr m_u_diff
Storage for F(d)-u_obs.
DesignVecGhostedPtr m_d
Current iterate of design parameter.
DesignVec & m_d0
A-priori estimate of design parameter.
std::shared_ptr< const Grid > m_grid
double m_val_design
Value of at the current iterate.
IPFunctional< array::Vector > & m_stateFunctional
Implementation of .
virtual std::shared_ptr< TerminationReason > formInitialGuess(Vec *v)
Callback from TaoBasicSolver to form the starting iterate for the minimization. See also.
virtual DesignVecPtr designSolution()
Value of , the solution of the minimization problem.
std::shared_ptr< typename ForwardProblem::DesignVecGhosted > DesignVecGhostedPtr
StateVec & m_u_obs
State parameter to match via F(d)=u_obs.
DesignVecPtr m_grad_design
Gradient of at the current iterate.
DesignVec m_dGlobal
Initial iterate of design parameter, stored without ghosts for the benefit of TAO.
StateVec m_adjointRHS
Temporary storage used in gradient computation.
ForwardProblem::DesignVecGhosted DesignVecGhosted
virtual void monitorTao(Tao tao)
Callback from TAO after each iteration. The call is forwarded to each element of our list of listener...
IPTaoTikhonovProblem(ForwardProblem &forward, DesignVec &d0, StateVec &u_obs, double eta, IPFunctional< DesignVec > &designFunctional, IPFunctional< StateVec > &stateFunctional)
virtual StateVecPtr stateSolution()
Final value of , where is the solution of the minimization.
Defines a Tikhonov minimization problem to be solved with a TaoBasicSolver.
Adaptor to connect a TAO objective function callback to a C++ object method.
Definition: TaoUtil.hh:363
Adaptor to connect a TAO monitoring callback to a C++ object method.
Definition: TaoUtil.hh:230
Adaptor to connect a TAO objective and gradient function callback to a C++ object method.
Definition: TaoUtil.hh:407
#define PISM_CHK(errcode, name)
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
static const double k
Definition: exactTestP.cc:42