tIP_SSATaucTaoTikhonovProblemLCL.hh - 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.hh (5951B)
---
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 #ifndef IP_SSATAUCTIKHONOVLCL_HH_39UGM4S2
20 #define IP_SSATAUCTIKHONOVLCL_HH_39UGM4S2
21
22 #include <memory>
23
24 #include <petscsys.h>
25
26 #include "TaoUtil.hh"
27 #include "IPTwoBlockVec.hh"
28 #include "pism/util/iceModelVec.hh"
29 #include "IP_SSATaucForwardProblem.hh"
30 #include "functional/IPFunctional.hh"
31
32 namespace pism {
33 namespace inverse {
34
35 class IP_SSATaucTaoTikhonovProblemLCL;
36
37 //! Iteration callback class for IP_SSATaucTaoTikhonovProblemLCL
38 /*! A class for objects receiving iteration callbacks from a
39 IP_SSATaucTaoTikhonovProblemLCL. These callbacks can be used to
40 monitor the solution, plot iterations, print diagnostic messages,
41 etc. IP_SSATaucTaoTikhonovProblemLCLListeners are ususally used via
42 a reference counted pointer
43 IP_SSATaucTaoTikhonovProblemLCLListeners::Ptr to allow for good
44 memory management when Listeners are created as subclasses of Python
45 classes.*/
46 class IP_SSATaucTaoTikhonovProblemLCLListener {
47 public:
48
49 typedef std::shared_ptr<IP_SSATaucTaoTikhonovProblemLCLListener> Ptr;
50
51 typedef IceModelVec2S DesignVec;
52 typedef IceModelVec2V StateVec;
53
54 IP_SSATaucTaoTikhonovProblemLCLListener() {}
55 virtual ~IP_SSATaucTaoTikhonovProblemLCLListener() {}
56
57 //! Callback called after each iteration.
58 //
59 // @param problem, The class calling the callback.
60 // @param eta Tikhonov penalty parameter.
61 // @param iter Current iteration count.
62 // @param objectiveValue Value of the state functional.
63 // @param designValue Value of the design functional.
64 // @param &d Value of the design variable.
65 // @param &diff_d Diference between design variable and a priori estimate.
66 // @param &grad_d Gradient of design functional
67 // @param &u Value of state variable
68 // @param &diff_u Difference between state variable and desired value.
69 // @param &grad_u Gradient of state functional
70 // @param constraints Residual for state variable being a solution of the %SSA
71 virtual void iteration(IP_SSATaucTaoTikhonovProblemLCL &problem,
72 double eta,
73 int iter,
74 double objectiveValue,
75 double designValue,
76 const DesignVec::Ptr &d,
77 const DesignVec::Ptr &diff_d,
78 const DesignVec::Ptr &grad_d,
79 const StateVec::Ptr &u,
80 const StateVec::Ptr &diff_u,
81 const StateVec::Ptr &grad_u,
82 const StateVec::Ptr &constraints) = 0;
83 };
84
85 //! \brief Defines a Tikhonov minimization problem of determining \f$\tau_c\f$ from %SSA velocities to be solved with a TaoBasicSolver using the tao_lcl algorithm.
86 /*! Experimental and not particularly functional. */
87 class IP_SSATaucTaoTikhonovProblemLCL {
88 public:
89 typedef IceModelVec2S DesignVec;
90 typedef IceModelVec2V StateVec;
91
92 typedef IP_SSATaucTaoTikhonovProblemLCLListener Listener;
93
94 IP_SSATaucTaoTikhonovProblemLCL(IP_SSATaucForwardProblem &ssaforward, DesignVec &d0, StateVec &u_obs, double eta,
95 IPFunctional<DesignVec> &designFunctional, IPFunctional<StateVec> &stateFunctional);
96
97 virtual ~IP_SSATaucTaoTikhonovProblemLCL();
98
99 virtual void addListener(Listener::Ptr listener) {
100 m_listeners.push_back(listener);
101 }
102
103 virtual StateVec::Ptr stateSolution();
104 virtual DesignVec::Ptr designSolution();
105
106 virtual void setInitialGuess(DesignVec &d0);
107
108 void connect(Tao tao);
109
110 void monitorTao(Tao tao);
111
112 virtual void evaluateObjectiveAndGradient(Tao tao, Vec x, double *value, Vec gradient);
113
114 virtual TerminationReason::Ptr formInitialGuess(Vec *x);
115
116 virtual void evaluateConstraints(Tao tao, Vec x, Vec r);
117
118 virtual void evaluateConstraintsJacobianState(Tao tao, Vec x, Mat Jstate, Mat Jpc, Mat Jinv,
119 MatStructure *s);
120
121 virtual void evaluateConstraintsJacobianDesign(Tao tao, Vec x, Mat Jdesign);
122
123 virtual void applyConstraintsJacobianDesign(Vec x, Vec y);
124 virtual void applyConstraintsJacobianDesignTranspose(Vec x, Vec y);
125
126 protected:
127
128 IP_SSATaucForwardProblem &m_ssaforward;
129
130 std::unique_ptr<IPTwoBlockVec> m_x;
131
132 DesignVec m_dGlobal;
133 DesignVec::Ptr m_d;
134 DesignVec &m_d0;
135 DesignVec::Ptr m_d_diff;
136 DesignVec m_dzeta;
137
138 StateVec::Ptr m_uGlobal;
139 StateVec m_u;
140 StateVec m_du;
141 StateVec &m_u_obs;
142 StateVec::Ptr m_u_diff;
143
144 DesignVec::Ptr m_grad_design;
145 StateVec::Ptr m_grad_state;
146
147 double m_eta;
148
149 double m_val_design;
150 double m_val_state;
151
152 StateVec::Ptr m_constraints;
153 petsc::Mat m_Jstate;
154 petsc::Mat m_Jdesign;
155
156 IceModelVec2S m_d_Jdesign;
157 IceModelVec2V m_u_Jdesign;
158
159 double m_constraintsScale;
160 double m_velocityScale;
161
162 IPFunctional<IceModelVec2S> &m_designFunctional;
163 IPFunctional<IceModelVec2V> &m_stateFunctional;
164
165 std::vector<Listener::Ptr> m_listeners;
166
167 static PetscErrorCode jacobian_design_callback(Mat A, Vec x, Vec y);
168 static PetscErrorCode jacobian_design_transpose_callback(Mat A, Vec x, Vec y);
169 };
170
171 } // end of namespace inverse
172 } // end of namespace pism
173
174 #endif /* end of include guard: IP_SSATAUCTIKHONOVLCL_HH_39UGM4S2 */