tIPDesignVariableParameterization.cc - 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
---
tIPDesignVariableParameterization.cc (6655B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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 #include <cmath>
20
21 #include "pism/util/iceModelVec.hh"
22 #include "IPDesignVariableParameterization.hh"
23 #include "pism/util/pism_options.hh"
24 #include "pism/util/ConfigInterface.hh"
25 #include "pism/util/IceGrid.hh"
26 #include "pism/util/error_handling.hh"
27
28 namespace pism {
29 namespace inverse {
30
31 //! Initializes the scale parameters of the parameterization.
32 /*! Every IPDesignVariableParameterization has an associated scale for the design variable
33 \f$d_{\rm scale}\f$ that equals 1 in internal units. The scale for a design variable named \a foo
34 is stored in an Config file as design_param_foo_scale. Subclasses may have additional
35 parameters that are follow the naming convention \a design_param_foo_*.
36
37 \param config The config file to read the scale parameters from.
38 \param design_var_name The associated name of the design variable, e.g. 'tauc' or 'hardav'
39 */
40 void IPDesignVariableParameterization::set_scales(const Config & config,
41 const std::string &design_var_name) {
42 std::string key("inverse.design.param_");
43 key += design_var_name;
44 key += "_scale";
45 m_d_scale = config.get_number(key);
46 }
47
48 //! Transforms a vector of \f$\zeta\f$ values to a vector of \f$d\f$ values.
49 void IPDesignVariableParameterization::convertToDesignVariable(IceModelVec2S &zeta,
50 IceModelVec2S &d,
51 bool communicate) {
52 PetscErrorCode ierr;
53
54 IceModelVec::AccessList list{&zeta, &d};
55
56 const IceGrid &grid = *zeta.grid();
57
58 ParallelSection loop(grid.com);
59 try {
60 for (Points p(grid); p; p.next()) {
61 const int i = p.i(), j = p.j();
62
63 this->toDesignVariable(zeta(i, j), &d(i, j), NULL);
64 if (std::isnan(d(i, j))) {
65 ierr = PetscPrintf(PETSC_COMM_WORLD,
66 "made a d nan zeta = %g d = %g\n",
67 zeta(i, j), d(i, j));
68 PISM_CHK(ierr, "PetscPrintf");
69 }
70 }
71 } catch (...) {
72 loop.failed();
73 }
74 loop.check();
75
76 if (communicate) {
77 d.update_ghosts();
78 }
79 }
80
81 //! Transforms a vector of \f$d\f$ values to a vector of \f$\zeta\f$ values.
82 void IPDesignVariableParameterization::convertFromDesignVariable(IceModelVec2S &d,
83 IceModelVec2S &zeta,
84 bool communicate) {
85 PetscErrorCode ierr;
86 IceModelVec::AccessList list{&zeta, &d};
87
88 const IceGrid &grid = *zeta.grid();
89
90 ParallelSection loop(grid.com);
91 try {
92 for (Points p(grid); p; p.next()) {
93 const int i = p.i(), j = p.j();
94
95 this->fromDesignVariable(d(i, j), &zeta(i, j));
96 if (std::isnan(zeta(i, j))) {
97 ierr = PetscPrintf(PETSC_COMM_WORLD,
98 "made a zeta nan d = %g zeta = %g\n",
99 d(i, j), zeta(i, j));
100 PISM_CHK(ierr, "PetscPrintf");
101 }
102 }
103 } catch (...) {
104 loop.failed();
105 }
106 loop.check();
107
108 if (communicate) {
109 zeta.update_ghosts();
110 }
111 }
112
113 void IPDesignVariableParamIdent::toDesignVariable(double p, double *value,
114 double *derivative) {
115 if (value != NULL) {
116 *value = m_d_scale*p;
117 }
118 if (derivative != NULL) {
119 *derivative = m_d_scale;
120 }
121 }
122
123 void IPDesignVariableParamIdent::fromDesignVariable(double d, double *OUTPUT) {
124 *OUTPUT = d / m_d_scale;
125 }
126
127
128 void IPDesignVariableParamSquare::toDesignVariable(double p, double *value,
129 double *derivative) {
130 if (value != NULL) {
131 *value = m_d_scale*p*p;
132 }
133 if (derivative != NULL) {
134 *derivative = m_d_scale*2*p;
135 }
136 }
137
138 void IPDesignVariableParamSquare::fromDesignVariable(double d, double *OUTPUT) {
139 if (d < 0) {
140 d = 0;
141 }
142 *OUTPUT = sqrt(d / m_d_scale);
143 }
144
145 void IPDesignVariableParamExp::set_scales(const Config &config, const std::string &design_var_name) {
146 IPDesignVariableParameterization::set_scales(config, design_var_name);
147
148 std::string key("inverse.design.param_");
149 key += design_var_name;
150 key += "_eps";
151 m_d_eps = config.get_number(key);
152 }
153
154 void IPDesignVariableParamExp::toDesignVariable(double p, double *value,
155 double *derivative) {
156 if (value != NULL) {
157 *value = m_d_scale*exp(p);
158 }
159
160 if (derivative != NULL) {
161 *derivative = m_d_scale*exp(p);
162 }
163 }
164
165 void IPDesignVariableParamExp::fromDesignVariable(double d, double *OUTPUT) {
166 if (d < m_d_eps) {
167 d = m_d_eps;
168 }
169 *OUTPUT = log(d / m_d_scale);
170 }
171
172
173 void IPDesignVariableParamTruncatedIdent::set_scales(const Config &config,
174 const std::string &design_var_name) {
175 IPDesignVariableParameterization::set_scales(config, design_var_name);
176
177 std::string key("inverse.design.param_trunc_");
178 key += design_var_name;
179 key += "0";
180
181 double d0 = config.get_number(key);
182 m_d0_sq = d0*d0 / (m_d_scale*m_d_scale);
183
184
185 key = "inverse.design.param_";
186 key += design_var_name;
187 key += "_eps";
188 m_d_eps = config.get_number(key);
189 }
190
191 void IPDesignVariableParamTruncatedIdent::toDesignVariable(double p,
192 double *value,
193 double *derivative) {
194 double alpha = sqrt(p*p + 4*m_d0_sq);
195 if (value != NULL) {
196 *value = m_d_scale*(p + alpha)*0.5;
197 }
198 if (derivative != NULL) {
199 *derivative = m_d_scale*(1 + p / alpha)*0.5;
200 }
201 }
202
203 void IPDesignVariableParamTruncatedIdent::fromDesignVariable(double d, double *OUTPUT) {
204 if (d < m_d_eps) {
205 d = m_d_eps;
206 }
207
208 double d_dimensionless = d / m_d_scale;
209 *OUTPUT = d_dimensionless - m_d0_sq / d_dimensionless;
210 }
211
212 } // end of namespace inverse
213 } // end of namespace pism