tIP_L2NormFunctional.cc - 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_L2NormFunctional.cc (8705B)
---
1 // Copyright (C) 2012, 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 #include "IP_L2NormFunctional.hh"
20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/pism_utilities.hh"
22
23 namespace pism {
24 namespace inverse {
25
26 void IP_L2NormFunctional2S::valueAt(IceModelVec2S &x, double *OUTPUT) {
27
28 const unsigned int Nq = m_quadrature.n();
29 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
30
31 // The value of the objective
32 double value = 0;
33
34 double x_q[Nq_max];
35
36 IceModelVec::AccessList list(x);
37
38 // Jacobian times weights for quadrature.
39 const double* W = m_quadrature.weights();
40
41 // Loop through all LOCAL elements.
42 const int
43 xs = m_element_index.lxs,
44 xm = m_element_index.lxm,
45 ys = m_element_index.lys,
46 ym = m_element_index.lym;
47
48 for (int j = ys; j < ys + ym; j++) {
49 for (int i = xs; i < xs + xm; i++) {
50 m_element.reset(i, j);
51
52 // Obtain values of x at the quadrature points for the element.
53 double tmp[fem::q1::n_chi];
54 m_element.nodal_values(x, tmp);
55 quadrature_point_values(m_quadrature, tmp, x_q);
56
57 for (unsigned int q = 0; q < Nq; q++) {
58 const double x_qq = x_q[q];
59 value += W[q]*x_qq*x_qq;
60 } // q
61 } // j
62 } // i
63
64 GlobalSum(m_grid->com, &value, OUTPUT, 1);
65 }
66
67 void IP_L2NormFunctional2S::dot(IceModelVec2S &a, IceModelVec2S &b, double *OUTPUT) {
68
69 const unsigned int Nq = m_quadrature.n();
70 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
71
72 // The value of the objective
73 double value = 0;
74
75 double a_q[Nq_max];
76 double b_q[Nq_max];
77
78 IceModelVec::AccessList list{&a, &b};
79
80 // Jacobian times weights for quadrature.
81 const double* W = m_quadrature.weights();
82
83 // Loop through all LOCAL elements.
84 const int
85 xs = m_element_index.lxs,
86 xm = m_element_index.lxm,
87 ys = m_element_index.lys,
88 ym = m_element_index.lym;
89
90 for (int j = ys; j < ys + ym; j++) {
91 for (int i = xs; i < xs + xm; i++) {
92 m_element.reset(i, j);
93
94 double tmp[fem::q1::n_chi];
95 m_element.nodal_values(a, tmp);
96 quadrature_point_values(m_quadrature, tmp, a_q);
97
98 m_element.nodal_values(b, tmp);
99 quadrature_point_values(m_quadrature, tmp, b_q);
100
101 for (unsigned int q = 0; q < Nq; q++) {
102 value += W[q]*a_q[q]*b_q[q];
103 } // q
104 } // j
105 } // i
106
107 GlobalSum(m_grid->com, &value, OUTPUT, 1);
108 }
109
110 void IP_L2NormFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S &gradient) {
111
112 const unsigned int Nk = fem::q1::n_chi;
113 const unsigned int Nq = m_quadrature.n();
114 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
115
116 // Clear the gradient before doing anything with it!
117 gradient.set(0);
118
119 double x_q[Nq_max];
120 double gradient_e[Nk];
121
122 IceModelVec::AccessList list{&x, &gradient};
123
124 // An Nq by Nk array of test function values.
125 const fem::Germs *test = m_quadrature.test_function_values();
126
127 // Jacobian times weights for quadrature.
128 const double* W = m_quadrature.weights();
129
130 // Loop through all local and ghosted elements.
131 const int
132 xs = m_element_index.xs,
133 xm = m_element_index.xm,
134 ys = m_element_index.ys,
135 ym = m_element_index.ym;
136
137 for (int j = ys; j < ys + ym; j++) {
138 for (int i = xs; i < xs + xm; i++) {
139
140 // Reset the DOF map for this element.
141 m_element.reset(i, j);
142
143 // Obtain values of x at the quadrature points for the element.
144 double tmp[Nk];
145 m_element.nodal_values(x, tmp);
146 quadrature_point_values(m_quadrature, tmp, x_q);
147
148 // Zero out the element-local residual in prep for updating it.
149 for (unsigned int k = 0; k < Nk; k++) {
150 gradient_e[k] = 0;
151 }
152
153 for (unsigned int q = 0; q < Nq; q++) {
154 const double x_qq = x_q[q];
155 for (unsigned int k = 0; k < Nk; k++) {
156 gradient_e[k] += 2*W[q]*x_qq*test[q][k].val;
157 } // k
158 } // q
159 m_element.add_contribution(gradient_e, gradient);
160 } // j
161 } // i
162 }
163
164 void IP_L2NormFunctional2V::valueAt(IceModelVec2V &x, double *OUTPUT) {
165
166 const unsigned int Nq = m_quadrature.n();
167 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
168
169 // The value of the objective
170 double value = 0;
171
172 Vector2 x_q[Nq_max];
173
174 IceModelVec::AccessList list(x);
175
176 // Jacobian times weights for quadrature.
177 const double* W = m_quadrature.weights();
178
179 // Loop through all local and ghosted elements.
180 const int
181 xs = m_element_index.lxs,
182 xm = m_element_index.lxm,
183 ys = m_element_index.lys,
184 ym = m_element_index.lym;
185
186 for (int j = ys; j < ys + ym; j++) {
187 for (int i = xs; i < xs + xm; i++) {
188 m_element.reset(i, j);
189
190 // Obtain values of x at the quadrature points for the element.
191 Vector2 tmp[fem::q1::n_chi];
192 m_element.nodal_values(x, tmp);
193 quadrature_point_values(m_quadrature, tmp, x_q);
194
195 for (unsigned int q = 0; q < Nq; q++) {
196 const Vector2 &x_qq = x_q[q];
197 value += W[q]*(x_qq.u*x_qq.u + x_qq.v*x_qq.v);
198 } // q
199 } // j
200 } // i
201
202 GlobalSum(m_grid->com, &value, OUTPUT, 1);
203 }
204
205 void IP_L2NormFunctional2V::dot(IceModelVec2V &a, IceModelVec2V &b, double *OUTPUT) {
206
207 const unsigned int Nq = m_quadrature.n();
208 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
209
210 // The value of the objective
211 double value = 0;
212
213 Vector2 a_q[Nq_max];
214 Vector2 b_q[Nq_max];
215
216 IceModelVec::AccessList list{&a, &b};
217
218 // Jacobian times weights for quadrature.
219 const double* W = m_quadrature.weights();
220
221 // Loop through all LOCAL elements.
222 const int
223 xs = m_element_index.lxs,
224 xm = m_element_index.lxm,
225 ys = m_element_index.lys,
226 ym = m_element_index.lym;
227
228 for (int j = ys; j < ys + ym; j++) {
229 for (int i = xs; i < xs + xm; i++) {
230 m_element.reset(i, j);
231
232 // Obtain values of x at the quadrature points for the element.
233 Vector2 tmp[fem::q1::n_chi];
234 m_element.nodal_values(a, tmp);
235 quadrature_point_values(m_quadrature, tmp, a_q);
236 m_element.nodal_values(b, tmp);
237 quadrature_point_values(m_quadrature, tmp, b_q);
238
239 for (unsigned int q = 0; q < Nq; q++) {
240 value += W[q]*(a_q[q].u*b_q[q].u + a_q[q].v*b_q[q].v);
241 } // q
242 } // j
243 } // i
244
245 GlobalSum(m_grid->com, &value, OUTPUT, 1);
246 }
247
248 void IP_L2NormFunctional2V::gradientAt(IceModelVec2V &x, IceModelVec2V &gradient) {
249
250 const unsigned int Nk = fem::q1::n_chi;
251 const unsigned int Nq = m_quadrature.n();
252 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
253
254 // Clear the gradient before doing anything with it!
255 gradient.set(0);
256
257 Vector2 x_q[Nq_max];
258 Vector2 gradient_e[Nk];
259
260 IceModelVec::AccessList list{&x, &gradient};
261
262 // An Nq by Nk array of test function values.
263 const fem::Germs *test = m_quadrature.test_function_values();
264
265 // Jacobian times weights for quadrature.
266 const double* W = m_quadrature.weights();
267
268 // Loop through all local and ghosted elements.
269 const int
270 xs = m_element_index.xs,
271 xm = m_element_index.xm,
272 ys = m_element_index.ys,
273 ym = m_element_index.ym;
274
275 for (int j = ys; j < ys + ym; j++) {
276 for (int i = xs; i < xs + xm; i++) {
277
278 // Reset the DOF map for this element.
279 m_element.reset(i, j);
280
281 // Obtain values of x at the quadrature points for the element.
282 Vector2 tmp[Nk];
283 m_element.nodal_values(x, tmp);
284 quadrature_point_values(m_quadrature, tmp, x_q);
285
286 // Zero out the element-local residual in prep for updating it.
287 for (unsigned int k = 0; k < Nk; k++) {
288 gradient_e[k].u = 0;
289 gradient_e[k].v = 0;
290 }
291
292 for (unsigned int q = 0; q < Nq; q++) {
293 const Vector2 &x_qq = x_q[q];
294 for (unsigned int k = 0; k < Nk; k++) {
295 double gcommon = 2*W[q]*test[q][k].val;
296 gradient_e[k].u += gcommon*x_qq.u;
297 gradient_e[k].v += gcommon*x_qq.v;
298 } // k
299 } // q
300 m_element.add_contribution(gradient_e, gradient);
301 } // j
302 } // i
303 }
304
305 } // end of namespace inverse
306 } // end of namespace pism