tIP_SSATaucForwardProblem.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_SSATaucForwardProblem.hh (8914B)
---
1 // Copyright (C) 2012, 2013, 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_SSATAUCFORWARDPROBLEM_HH_4AEVR4Z
20 #define IP_SSATAUCFORWARDPROBLEM_HH_4AEVR4Z
21
22 #include "pism/stressbalance/ssa/SSAFEM.hh"
23 #include "IPDesignVariableParameterization.hh"
24 #include "pism/util/petscwrappers/KSP.hh"
25 #include "pism/util/petscwrappers/Mat.hh"
26
27 namespace pism {
28 namespace inverse {
29
30 //! Implements the forward problem of the map taking \f$\tau_c\f$ to the corresponding solution of the %SSA.
31 /*! The class SSAFEM solves the %SSA, and the solution depends on a large number of parameters. Considering
32 all of these to be fixed except for \f$\tau_c\f$, we obtain a map \f$F_{\rm SSA}\f$ taking
33 \f$\tau_c\f$ to the corresponding solution \f$u_{\rm SSA}\f$ of the %SSA. This is a forward problem which
34 we would like to invert; given \f$u_{\rm SSA}\f$ determine \f$\tau_c\f$ such that
35 \f$F_{\rm SSA}(\tau_c) = u_{\rm SSA}\f$. The forward problem actually implemented by
36 IP_SSATaucForwardProblem is a mild variation \f$F_{\rm SSA}\f$.
37
38 First, given the constraint \f$\tau_c\ge 0\f$ it can be helpful to parameterize \f$\tau_c\f$ by some other
39 parameter \f$\zeta\f$,
40 \f[
41 \tau_c = g(\zeta),
42 \f]
43 where the function \f$g\f$ is non-negative. The function \f$g\f$ is specified by an instance
44 of IPDesignVariableParameterization.
45
46 Second, there may be locations where the value of \f$\tau_c\f$ (and hence \f$\zeta\f$)
47 is known a-priori, and should not be adjusted. Let \f$\pi\f$ be the map that replaces
48 the values of \f$\zeta\f$ with known values at these locations.
49
50 IP_SSATaucForwardProblem implements the forward problem
51 \f[
52 F(\zeta) = F_{\rm SSA}(g(\pi(\zeta))).
53 \f]
54
55 Algorithms to solve the inverse problem make use of variations on the linearization
56 of \f$F\f$, which are explained below.
57
58 The solution of the %SSA in SSAFEM is implemented using SNES to solve
59 a nonlinear residual function of the form
60 \f[
61 \mathcal{R}_{SSA}(u;\tau_c, \text{other parameters})= \vec 0,
62 \f]
63 usually using Newton's method. Let
64 \f[
65 \mathcal{R}(u, \zeta) = \mathcal{R}_{SSA}(u; g(\pi(\zeta)), \text{other parameters}).
66 \f]
67
68 The derivative of \f$\mathcal{R}\f$ with respect to \f$ u\f$ is called the state Jacobian,
69 \f$J_{\rm State}\f$. Specifically, if \f$u=[U_j]\f$ then
70 \f[
71 (J_{\rm State})_{ij} = \frac{d \mathcal{R}_i}{dU_j}.
72 \f]
73 This is exactly the same Jacobian that is computed by SSAFEM for solving the %SSA via SNES. The
74 matrix for the design Jacobian can be obtained with \ref assemble_jacobian_state.
75
76 The derivative of \f$\mathcal{R}\f$ with respect to \f$ \zeta\f$ is called the design Jacobian,
77 \f$J_{\rm Design}\f$. Specifically, if \f$\zeta=[Z_k]\f$ then
78 \f[
79 (J_{\rm Design})_{ik} = \frac{d \mathcal{R}_i}{dZ_k}.
80 \f]
81 The map \f$J_{\rm Design}\f$ can be applied to a vector \f$d\zeta\f$ with
82 apply_jacobian_design. For inverse methods using adjoints, one also
83 needs to be able to apply the transpose of \f$J_{\rm Design}\f$,
84 which is done using \ref apply_jacobian_design_transpose.
85
86 The forward problem map \f$F\f$ solves the implicit equation
87 \f[
88 \mathcal{R}(F(\zeta), \zeta) = 0.
89 \f]
90 Its linearisation \f$DF\f$ then satisfies
91 \f[
92 J_{\rm State}\; DF\; d\zeta + J_{\rm Design}\; d\zeta = 0
93 \f]
94 for any perturbation \f$d\zeta\f$. Hence
95 \f[
96 DF = -J_{\rm State}^{-1}\circ J_{\rm Design}.
97 \f]
98 This derivative is sometimes called the reduced gradient in the literature. To
99 apply \f$DF\f$ to a perturbation \f$d\zeta\f$, use \ref apply_linearization. Adjoint
100 methods require the transpose of this map; to apply \f$DF^t\f$ to \f$du\f$ use
101 \ref apply_linearization_transpose.
102 */
103 class IP_SSATaucForwardProblem : public stressbalance::SSAFEM
104 {
105 public:
106
107 typedef IceModelVec2S DesignVec; ///< The function space for the design variable, i.e. \f$\tau_c\f$.
108 typedef IceModelVec2V StateVec; ///< The function space for the state variable, \f$u_{\rm SSA}\f$.
109
110 //! Constructs from the same objects as SSAFEM, plus a specification of how \f$\tau_c\f$ is parameterized.
111 IP_SSATaucForwardProblem(IceGrid::ConstPtr g,
112 IPDesignVariableParameterization &tp);
113
114 virtual ~IP_SSATaucForwardProblem();
115
116 void init();
117
118 //! Selects nodes where \f$\tau_c\f$ (more specifically \f$\zeta\f$) should not be adjusted.
119 /*! The paramter \a locations should be set to 1 at each node where \f$\tau_c\f$
120 is fixed. The forward map then effectively treats the design space as the subspace
121 of nodes where \a locations is 0. Tangent vectors to this subspace, as would be
122 generated by, e.g., \f$J_{\rm Design}^t\f$ are represented as vectors in the full
123 space with entries set to zero in the fixed locations. These can safely be added
124 to preexisting values of \f$\zeta\f$ without changing the entries of \f$\zeta\f$ at the
125 fixed locations. Inversion can be done by setting an initial value of \f$\zeta\f$
126 having the desired values in the fixed locations, and using set_tauc_fixed_locations()
127 to indicate the nodes that should not be changed.
128 */
129 virtual void set_tauc_fixed_locations(IceModelVec2Int &locations)
130 {
131 m_fixed_tauc_locations = &locations;
132 }
133
134 //! Returns the last solution of the %SSA as computed by \ref linearize_at.
135 virtual IceModelVec2V::Ptr solution() {
136 m_velocity_shared->copy_from(m_velocity);
137 return m_velocity_shared;
138 }
139
140 //! Exposes the \f$\tau_c\f$ parameterization being used.
141 virtual IPDesignVariableParameterization & tauc_param() {
142 return m_tauc_param;
143 }
144
145 virtual void set_design(IceModelVec2S &zeta);
146
147 virtual TerminationReason::Ptr linearize_at(IceModelVec2S &zeta);
148
149 virtual void assemble_residual(IceModelVec2V &u, IceModelVec2V &R);
150 virtual void assemble_residual(IceModelVec2V &u, Vec R);
151
152 virtual void assemble_jacobian_state(IceModelVec2V &u, Mat J);
153
154 virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, IceModelVec2V &du);
155 virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, Vec du);
156 virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, Vector2 **du_a);
157
158 virtual void apply_jacobian_design_transpose(IceModelVec2V &u, IceModelVec2V &du, IceModelVec2S &dzeta);
159 virtual void apply_jacobian_design_transpose(IceModelVec2V &u, IceModelVec2V &du, Vec dzeta);
160 virtual void apply_jacobian_design_transpose(IceModelVec2V &u, IceModelVec2V &du, double **dzeta);
161
162 virtual void apply_linearization(IceModelVec2S &dzeta, IceModelVec2V &du);
163 virtual void apply_linearization_transpose(IceModelVec2V &du, IceModelVec2S &dzeta);
164
165 //! Exposes the DMDA of the underlying grid for the benefit of TAO.
166 virtual void get_da(DM *da) {
167 *da = *m_da;
168 }
169
170 protected:
171
172 /// Current value of zeta, provided from caller.
173 IceModelVec2S *m_zeta;
174 /// Storage for d_zeta with ghosts, if needed when an argument d_zeta is ghost-less.
175 IceModelVec2S m_dzeta_local;
176 /// Storage for tauc (avoids modifying fields obtained via pism::Vars)
177 IceModelVec2S m_tauc_copy;
178
179 /// Locations where \f$\tau_c\f$ should not be adjusted.
180 IceModelVec2Int *m_fixed_tauc_locations;
181
182 /// The function taking \f$\zeta\f$ to \f$\tau_c\f$.
183 IPDesignVariableParameterization &m_tauc_param;
184
185 /// Copy of the velocity field managed using a shared pointer.
186 IceModelVec2V::Ptr m_velocity_shared;
187
188 /// Temporary storage when state vectors need to be used without ghosts.
189 IceModelVec2V m_du_global;
190 /// Temporary storage when state vectors need to be used with ghosts.
191 IceModelVec2V m_du_local;
192
193 fem::ElementIterator m_element_index;
194 fem::ElementMap m_element;
195 fem::Q1Quadrature4 m_quadrature;
196
197 /// KSP used in \ref apply_linearization and \ref apply_linearization_transpose
198 petsc::KSP m_ksp;
199 /// Mat used in \ref apply_linearization and \ref apply_linearization_transpose
200 petsc::Mat m_J_state;
201
202 SNESConvergedReason m_reason;
203
204 /// Flag indicating that the state jacobian matrix needs rebuilding.
205 bool m_rebuild_J_state;
206 };
207
208 } // end of namespace inverse
209 } // end of namespace pism
210
211 #endif /* end of include guard: IP_SSATAUCFORWARDPROBLEM_HH_4AEVR4Z */