tIP_SSAHardavForwardProblem.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_SSAHardavForwardProblem.hh (8918B)
---
1 // Copyright (C) 2013, 2014, 2015, 2016, 2017 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 2 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_SSAHARDAVFORWARDPROBLEM_HH_HAD68BK7
20 #define IP_SSAHARDAVFORWARDPROBLEM_HH_HAD68BK7
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_SSAHardavForwardProblem 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_SSAHardavForwardProblem 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_SSAHardavForwardProblem : 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_SSAHardavForwardProblem(IceGrid::ConstPtr g,
112 IPDesignVariableParameterization &tp);
113
114 virtual ~IP_SSAHardavForwardProblem();
115
116 //! Selects nodes where \f$\tau_c\f$ (more specifically \f$\zeta\f$) should not be adjusted.
117 /*! The paramter \a locations should be set to 1 at each node where \f$\tau_c\f$
118 is fixed. The forward map then effectively treats the design space as the subspace
119 of nodes where \a locations is 0. Tangent vectors to this subspace, as would be
120 generated by, e.g., \f$J_{\rm Design}^t\f$ are represented as vectors in the full
121 space with entries set to zero in the fixed locations. These can safely be added
122 to preexisting values of \f$\zeta\f$ without changing the entries of \f$\zeta\f$ at the
123 fixed locations. Inversion can be done by setting an initial value of \f$\zeta\f$
124 having the desired values in the fixed locations, and using set_tauc_fixed_locations()
125 to indicated the nodes that should not be changed.
126 */
127 virtual void set_design_fixed_locations(IceModelVec2Int &locations)
128 {
129 m_fixed_design_locations = &locations;
130 }
131
132 //! Returns the last solution of the %SSA as computed by \ref linearize_at.
133 virtual IceModelVec2V::Ptr solution() {
134 m_velocity_shared->copy_from(m_velocity);
135 return m_velocity_shared;
136 }
137
138 //! Exposes the design variable parameterization being used.
139 virtual IPDesignVariableParameterization & design_param() {
140 return m_design_param;
141 }
142
143 void init();
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 IceModelVec2S *m_zeta; ///< Current value of zeta, provided from caller.
173 IceModelVec2S m_dzeta_local; ///< Storage for d_zeta with ghosts, if needed when an argument d_zeta is ghost-less.
174
175 IceModelVec2Int *m_fixed_design_locations; ///< Locations where \f$\tau_c\f$ should not be adjusted.
176
177 IPDesignVariableParameterization &m_design_param; ///< The function taking \f$\zeta\f$ to \f$\tau_c\f$.
178
179 IceModelVec2V::Ptr m_velocity_shared;
180
181 IceModelVec2V m_du_global; ///< Temporary storage when state vectors need to be used without ghosts.
182 IceModelVec2V m_du_local; ///< Temporary storage when state vectors need to be used with ghosts.
183 IceModelVec2S m_hardav;
184
185 fem::ElementIterator m_element_index;
186 fem::ElementMap m_element;
187 fem::Q1Quadrature4 m_quadrature;
188
189 petsc::KSP m_ksp; ///< KSP used in \ref apply_linearization and \ref apply_linearization_transpose
190 petsc::Mat m_J_state; ///< Mat used in \ref apply_linearization and \ref apply_linearization_transpose
191
192 SNESConvergedReason m_reason;
193
194 bool m_rebuild_J_state; ///< Flag indicating that the state jacobian matrix needs rebuilding.
195 };
196
197 } // end of namespace inverse
198 } // end of namespace pism
199
200 #endif /* end of include guard: IP_SSAHARDAVFORWARDPROBLEM_HH_HAD68BK7 */