tIP_SSATaucTikhonovGNSolver.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_SSATaucTikhonovGNSolver.hh (4795B)
---
1 // Copyright (C) 2012, 2014, 2015, 2016, 2017, 2019 David Maxwell
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_SSATAUCTIKHONOVGN_HH_SIU7F33G
20 #define IP_SSATAUCTIKHONOVGN_HH_SIU7F33G
21
22 #include "IP_SSATaucForwardProblem.hh"
23 #include "pism/util/TerminationReason.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/iceModelVec.hh"
26 #include "pism/util/petscwrappers/KSP.hh"
27 #include "functional/IPFunctional.hh"
28
29 namespace pism {
30 namespace inverse {
31
32 template<class C, void (C::*MultiplyCallback)(Vec,Vec) >
33 class MatrixMultiplyCallback {
34 public:
35 static void connect(Mat A) {
36 PetscErrorCode ierr;
37 ierr = MatShellSetOperation(A, MATOP_MULT,
38 (void(*)(void))MatrixMultiplyCallback::callback);
39 PISM_CHK(ierr, "MatShellSetOperation");
40 }
41 protected:
42 static PetscErrorCode callback(Mat A, Vec x, Vec y) {
43 try {
44 C *ctx;
45 PetscErrorCode ierr = MatShellGetContext(A,&ctx);
46 PISM_CHK(ierr, "MatShellGetContext");
47 (ctx->*MultiplyCallback)(x,y);
48 } catch (...) {
49 MPI_Comm com = MPI_COMM_SELF;
50 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
51 handle_fatal_errors(com);
52 SETERRQ(com, 1, "A PISM callback failed");
53 }
54 return 0;
55 }
56 };
57
58 class IP_SSATaucTikhonovGNSolver {
59 public:
60 typedef IceModelVec2S DesignVec;
61 typedef IceModelVec2V StateVec;
62 // typedef IP_SSATaucTikhonovGNSolverListener Listener;
63
64 IP_SSATaucTikhonovGNSolver(IP_SSATaucForwardProblem &ssaforward, DesignVec &d0, StateVec &u_obs, double eta,
65 IPInnerProductFunctional<DesignVec> &designFunctional, IPInnerProductFunctional<StateVec> &stateFunctional);
66
67 ~IP_SSATaucTikhonovGNSolver();
68
69 virtual StateVec::Ptr stateSolution() {
70 return m_ssaforward.solution();
71 }
72
73 virtual DesignVec::Ptr designSolution() {
74 return m_d;
75 }
76
77 virtual void setInitialGuess(DesignVec &d) {
78 m_d->copy_from(d);
79 }
80
81 //! Sets the desired target misfit (in units of \f$\sqrt{J_{\rm misfit}}\f$).
82 virtual void setTargetMisfit(double misfit) {
83 m_target_misfit = misfit;
84 }
85
86 virtual void evaluateGNFunctional(DesignVec &h, double *value);
87
88 virtual void apply_GN(IceModelVec2S &h, IceModelVec2S &out);
89 virtual void apply_GN(Vec h, Vec out);
90
91 virtual TerminationReason::Ptr init();
92
93 virtual TerminationReason::Ptr check_convergence();
94
95 virtual TerminationReason::Ptr solve();
96
97 virtual TerminationReason::Ptr evaluate_objective_and_gradient();
98
99 protected:
100
101 virtual void assemble_GN_rhs(DesignVec &out);
102
103 virtual TerminationReason::Ptr solve_linearized();
104
105 virtual TerminationReason::Ptr compute_dlogalpha(double *dalpha);
106
107 virtual TerminationReason::Ptr linesearch();
108
109 const unsigned int m_design_stencil_width;
110 const unsigned int m_state_stencil_width;
111
112 IP_SSATaucForwardProblem &m_ssaforward;
113
114 DesignVec m_x;
115 DesignVec m_y;
116
117 DesignVec m_tmp_D1Global;
118 DesignVec m_tmp_D2Global;
119 DesignVec m_tmp_D1Local;
120 DesignVec m_tmp_D2Local;
121 StateVec m_tmp_S1Global;
122 StateVec m_tmp_S2Global;
123 StateVec m_tmp_S1Local;
124 StateVec m_tmp_S2Local;
125
126 DesignVec m_GN_rhs;
127
128 DesignVec::Ptr m_d;
129 DesignVec &m_d0;
130 DesignVec m_dGlobal;
131 DesignVec m_d_diff;
132 DesignVec m_d_diff_lin;
133
134 DesignVec m_h;
135 DesignVec m_hGlobal;
136 DesignVec m_dalpha_rhs;
137 DesignVec m_dh_dalpha;
138 DesignVec m_dh_dalphaGlobal;
139
140 DesignVec m_grad_design;
141 DesignVec m_grad_state;
142 DesignVec m_gradient;
143
144 double m_val_design, m_val_state, m_value;
145
146 StateVec &m_u_obs;
147 StateVec m_u_diff;
148
149 petsc::KSP m_ksp;
150 petsc::Mat m_mat_GN;
151
152 double m_eta;
153 IPInnerProductFunctional<DesignVec> &m_designFunctional;
154 IPInnerProductFunctional<StateVec> &m_stateFunctional;
155
156 double m_alpha;
157 double m_logalpha;
158 double m_target_misfit;
159
160 int m_iter, m_iter_max;
161 bool m_tikhonov_adaptive;
162 double m_vel_scale;
163 double m_tikhonov_rtol, m_tikhonov_atol, m_tikhonov_ptol;
164
165 MPI_Comm m_comm;
166 Logger::ConstPtr m_log;
167 };
168
169 } // end of namespace inverse
170 } // end of namespace pism
171
172 #endif /* end of include guard: IP_SSATAUCTIKHONOVGN_HH_SIU7F33G */