tIPFunctional.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tIPFunctional.hh (4750B)
---
1 // Copyright (C) 2012, 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 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 IPFUNCTIONAL_HH_1E2DIXE6
20 #define IPFUNCTIONAL_HH_1E2DIXE6
21
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/FETools.hh"
24
25 namespace pism {
26
27 //! Inverse modeling code.
28 namespace inverse {
29 //! Abstract base class for functions from ice model vectors to \f$\mathbb{R}\f$.
30 /*! Inverse problems frequently involve minimizing a functional,
31 such such as the misfit
32 \f[
33 J(u) = \int_\Omega |u-u_{\rm obs}|^2
34 \f]
35 between observed (\f$u_{\rm obs}\f$) and modeled (\f$u\f$)
36 surface velocities. Subclasses of IPFunctional define such maps,
37 and permit computation of their gradients.
38 */
39 template<class IMVecType>
40 class IPFunctional {
41
42 public:
43 IPFunctional(IceGrid::ConstPtr grid)
44 : m_grid(grid),
45 m_element_index(*m_grid),
46 m_element(*m_grid),
47 m_quadrature(grid->dx(), grid->dy(), 1.0)
48 {
49 }
50
51 virtual ~IPFunctional() {};
52
53 //! Computes the value of the functional at the vector x.
54 virtual void valueAt(IMVecType &x, double *OUTPUT) = 0;
55
56 //! Computes the gradient of the functional at the vector x.
57 /*! On an \f$m\times n\f$ IceGrid, an IceModelVec \f$x\f$ with \f$d\f$
58 degrees of freedom will be \f$d m n\f$-dimensional with components \f$x_i\f$.
59 The gradient computed here is the vector of directional derivatives \f$\nabla J\f$ of the functional
60 \f$J\f$ with respect to \f$x\f$. Concretely, the \f$i^{\rm th}\f$ component of \f$\nabla J\f$
61 is
62 \f[
63 \nabla J_i = \frac{\partial}{\partial x_i} J(x).
64 \f]
65 This vector is returned as `gradient`.
66 */
67 virtual void gradientAt(IMVecType &x, IMVecType &gradient) = 0;
68
69 protected:
70 IceGrid::ConstPtr m_grid;
71
72 fem::ElementIterator m_element_index;
73 fem::ElementMap m_element;
74 fem::Q1Quadrature4 m_quadrature;
75
76 private:
77 // Hide copy/assignment operations
78 IPFunctional(IPFunctional const &);
79 IPFunctional & operator=(IPFunctional const &);
80
81 };
82
83 //! Abstract base class for IPFunctionals arising from an inner product.
84 /*!
85 Frequently functionals have the structure
86 \f[
87 J(u) = Q(u, u)
88 \f]
89 where \f$Q\f$ is a symmetric, non-negative definite, bilinear form. Certain
90 minimization algorithms only apply to such functionals, which should subclass
91 from IPInnerProductFunctional.
92 */
93 template<class IMVecType>
94 class IPInnerProductFunctional : public IPFunctional<IMVecType> {
95
96 public:
97 IPInnerProductFunctional(IceGrid::ConstPtr grid) : IPFunctional<IMVecType>(grid) {};
98
99 //! Computes the inner product \f$Q(a, b)\f$.
100 virtual void dot(IMVecType &a, IMVecType &b, double *OUTPUT) = 0;
101
102 //! Computes the interior product of a vector with the IPInnerProductFunctional's underlying bilinear form.
103 /*! If \f$Q(x, y)\f$ is a bilinear form, and \f$a\f$ is a vector, then the
104 interior product of \f$a\f$ with \f$Q\f$ is the functional
105 \f[
106 I(z) = Q(a, z).
107 \f]
108 Such a functional is always linear and hence can be represented by taking
109 the standard dot product with some vector \f$y\f$:
110 \f[
111 I(z) = y^T z.
112 \f]
113 This method returns the vector \f$y\f$.
114 */
115 virtual void interior_product(IMVecType &x, IMVecType &y) {
116 this->gradientAt(x, y);
117 y.scale(0.5);
118 }
119
120 //! Assembles the matrix \f$Q_{ij}\f$ corresponding to the bilinear form.
121 /*! If \f$\{x_i\}\f$ is a basis for the vector space IMVecType,
122 \f$Q_{ij}= Q(x_i, x_j)\f$. */
123 // virtual void assemble_form(Mat Q) = 0;
124
125 };
126
127 //! Computes finite difference approximations of a IPFunctional<IceModelVec2S> gradient.
128 /*! Useful for debugging a hand coded gradient. */
129 void gradientFD(IPFunctional<IceModelVec2S> &f, IceModelVec2S &x, IceModelVec2S &gradient);
130
131 //! Computes finite difference approximations of a IPFunctional<IceModelVec2V> gradient.
132 /*! Useful for debugging a hand coded gradient. */
133 void gradientFD(IPFunctional<IceModelVec2V> &f, IceModelVec2V &x, IceModelVec2V &gradient);
134
135 } // end of namespace inverse
136 } // end of namespace pism
137
138 #endif /* end of include guard: FUNCTIONAL_HH_1E2DIXE6 */