tIPGroundedIceH1NormFunctional.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
---
tIPGroundedIceH1NormFunctional.cc (9643B)
---
1 // Copyright (C) 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 #include "IPGroundedIceH1NormFunctional.hh"
20 #include "pism/util/error_handling.hh"
21 #include "pism/util/IceGrid.hh"
22 #include "pism/util/pism_utilities.hh"
23 #include "pism/util/IceModelVec2CellType.hh"
24
25 namespace pism {
26 namespace inverse {
27
28 void IPGroundedIceH1NormFunctional2S::valueAt(IceModelVec2S &x, double *OUTPUT) {
29
30 const unsigned int Nk = fem::q1::n_chi;
31 const unsigned int Nq = m_quadrature.n();
32 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
33
34 // The value of the objective
35 double value = 0;
36
37 double x_e[Nk];
38 double x_q[Nq_max], dxdx_q[Nq_max], dxdy_q[Nq_max];
39
40 IceModelVec::AccessList list{&x, &m_ice_mask};
41
42 // Jacobian times weights for quadrature.
43 const double* W = m_quadrature.weights();
44
45 fem::DirichletData_Scalar dirichletBC(m_dirichletIndices, NULL);
46
47 // Loop through all LOCAL elements.
48 const int
49 xs = m_element_index.lxs,
50 xm = m_element_index.lxm,
51 ys = m_element_index.lys,
52 ym = m_element_index.lym;
53
54 for (int j=ys; j<ys+ym; j++) {
55 for (int i=xs; i<xs+xm; i++) {
56 bool all_grounded_ice = (m_ice_mask.grounded_ice(i, j) and
57 m_ice_mask.grounded_ice(i+1, j) and
58 m_ice_mask.grounded_ice(i, j+1) and
59 m_ice_mask.grounded_ice(i+1, j+1));
60
61 if (! all_grounded_ice) {
62 continue;
63 }
64
65 m_element.reset(i, j);
66
67 // Obtain values of x at the quadrature points for the element.
68 m_element.nodal_values(x, x_e);
69 if (dirichletBC) {
70 dirichletBC.enforce_homogeneous(m_element, x_e);
71 }
72 quadrature_point_values(m_quadrature, x_e, x_q, dxdx_q, dxdy_q);
73
74 for (unsigned int q=0; q<Nq; q++) {
75 value += W[q]*(m_cL2*x_q[q]*x_q[q]+ m_cH1*(dxdx_q[q]*dxdx_q[q]+dxdy_q[q]*dxdy_q[q]));
76 } // q
77 } // j
78 } // i
79
80 GlobalSum(m_grid->com, &value, OUTPUT, 1);
81 }
82
83 void IPGroundedIceH1NormFunctional2S::dot(IceModelVec2S &a, IceModelVec2S &b, double *OUTPUT) {
84
85 const unsigned int Nk = fem::q1::n_chi;
86 const unsigned int Nq = m_quadrature.n();
87 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
88
89 // The value of the objective
90 double value = 0;
91
92 double a_e[Nk];
93 double a_q[Nq_max], dadx_q[Nq_max], dady_q[Nq_max];
94
95 IceModelVec::AccessList list{&a, &b, &m_ice_mask};
96
97 double b_e[Nk];
98 double b_q[Nq_max], dbdx_q[Nq_max], dbdy_q[Nq_max];
99
100 // Jacobian times weights for quadrature.
101 const double* W = m_quadrature.weights();
102
103 fem::DirichletData_Scalar dirichletBC(m_dirichletIndices, NULL);
104
105 // Loop through all LOCAL elements.
106 const int
107 xs = m_element_index.lxs,
108 xm = m_element_index.lxm,
109 ys = m_element_index.lys,
110 ym = m_element_index.lym;
111
112 for (int j=ys; j<ys+ym; j++) {
113 for (int i=xs; i<xs+xm; i++) {
114 bool all_grounded_ice = (m_ice_mask.grounded_ice(i, j) and
115 m_ice_mask.grounded_ice(i+1, j) and
116 m_ice_mask.grounded_ice(i, j+1) and
117 m_ice_mask.grounded_ice(i+1, j+1));
118
119 if (! all_grounded_ice) {
120 continue;
121 }
122
123 m_element.reset(i, j);
124
125 // Obtain values of x at the quadrature points for the element.
126 m_element.nodal_values(a, a_e);
127 if (dirichletBC) {
128 dirichletBC.enforce_homogeneous(m_element, a_e);
129 }
130 quadrature_point_values(m_quadrature, a_e, a_q, dadx_q, dady_q);
131
132 m_element.nodal_values(b, b_e);
133 if (dirichletBC) {
134 dirichletBC.enforce_homogeneous(m_element, b_e);
135 }
136 quadrature_point_values(m_quadrature, b_e, b_q, dbdx_q, dbdy_q);
137
138 for (unsigned int q=0; q<Nq; q++) {
139 value += W[q]*(m_cL2*a_q[q]*b_q[q]+ m_cH1*(dadx_q[q]*dbdx_q[q]+dady_q[q]*dbdy_q[q]));
140 } // q
141 } // j
142 } // i
143
144 GlobalSum(m_grid->com, &value, OUTPUT, 1);
145 }
146
147
148 void IPGroundedIceH1NormFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S &gradient) {
149
150 const unsigned int Nk = fem::q1::n_chi;
151 const unsigned int Nq = m_quadrature.n();
152 const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
153
154 // Clear the gradient before doing anything with it!
155 gradient.set(0);
156
157 double x_e[Nk];
158 double x_q[Nq_max], dxdx_q[Nq_max], dxdy_q[Nq_max];
159
160 IceModelVec::AccessList list{&x, &gradient, &m_ice_mask};
161
162 double gradient_e[Nk];
163
164 // An Nq by Nk array of test function values.
165 const fem::Germs *test = m_quadrature.test_function_values();
166
167 // Jacobian times weights for quadrature.
168 const double* W = m_quadrature.weights();
169
170 fem::DirichletData_Scalar dirichletBC(m_dirichletIndices, NULL);
171
172 // Loop through all local and ghosted elements.
173 const int
174 xs = m_element_index.xs,
175 xm = m_element_index.xm,
176 ys = m_element_index.ys,
177 ym = m_element_index.ym;
178
179 for (int j=ys; j<ys+ym; j++) {
180 for (int i=xs; i<xs+xm; i++) {
181 bool all_grounded_ice = (m_ice_mask.grounded_ice(i, j) and
182 m_ice_mask.grounded_ice(i+1, j) and
183 m_ice_mask.grounded_ice(i, j+1) and
184 m_ice_mask.grounded_ice(i+1, j+1));
185
186 if (! all_grounded_ice) {
187 continue;
188 }
189
190 // Reset the DOF map for this element.
191 m_element.reset(i, j);
192
193 // Obtain values of x at the quadrature points for the element.
194 m_element.nodal_values(x, x_e);
195 if (dirichletBC) {
196 dirichletBC.constrain(m_element);
197 dirichletBC.enforce_homogeneous(m_element, x_e);
198 }
199 quadrature_point_values(m_quadrature, x_e, x_q, dxdx_q, dxdy_q);
200
201 // Zero out the element-local residual in prep for updating it.
202 for (unsigned int k=0; k<Nk; k++) {
203 gradient_e[k] = 0;
204 }
205
206 for (unsigned int q=0; q<Nq; q++) {
207 const double &x_qq=x_q[q];
208 const double &dxdx_qq=dxdx_q[q], &dxdy_qq=dxdy_q[q];
209 for (unsigned int k=0; k<Nk; k++) {
210 gradient_e[k] += 2*W[q]*(m_cL2*x_qq*test[q][k].val +
211 m_cH1*(dxdx_qq*test[q][k].dx + dxdy_qq*test[q][k].dy));
212 } // k
213 } // q
214 m_element.add_contribution(gradient_e, gradient);
215 } // j
216 } // i
217 }
218
219 void IPGroundedIceH1NormFunctional2S::assemble_form(Mat form) {
220
221 const unsigned int Nk = fem::q1::n_chi;
222 const unsigned int Nq = m_quadrature.n();
223
224 PetscErrorCode ierr;
225 int i, j;
226
227 // Zero out the Jacobian in preparation for updating it.
228 ierr = MatZeroEntries(form);
229 PISM_CHK(ierr, "MatZeroEntries");
230
231 // Jacobian times weights for quadrature.
232 const double* W = m_quadrature.weights();
233
234 fem::DirichletData_Scalar zeroLocs(m_dirichletIndices, NULL);
235
236 IceModelVec::AccessList list(m_ice_mask);
237
238 // Values of the finite element test functions at the quadrature points.
239 // This is an Nq by Nk array of function germs (Nq=#of quad pts, Nk=#of test functions).
240 const fem::Germs *test = m_quadrature.test_function_values();
241
242 // Loop through all the elements.
243 const int
244 xs = m_element_index.xs,
245 xm = m_element_index.xm,
246 ys = m_element_index.ys,
247 ym = m_element_index.ym;
248
249 for (j=ys; j<ys+ym; j++) {
250 for (i=xs; i<xs+xm; i++) {
251 bool all_grounded_ice = (m_ice_mask.grounded_ice(i, j) and
252 m_ice_mask.grounded_ice(i+1, j) and
253 m_ice_mask.grounded_ice(i, j+1) and
254 m_ice_mask.grounded_ice(i+1, j+1));
255 if (! all_grounded_ice) {
256 continue;
257 }
258
259 // Element-local Jacobian matrix (there are Nk vector valued degrees
260 // of freedom per elment, for a total of (2*Nk)*(2*Nk) = 16
261 // entries in the local Jacobian.
262 double K[Nk][Nk];
263
264 // Initialize the map from global to local degrees of freedom for this element.
265 m_element.reset(i, j);
266
267 // Don't update rows/cols where we project to zero.
268 if (zeroLocs) {
269 zeroLocs.constrain(m_element);
270 }
271
272 // Build the element-local Jacobian.
273 ierr = PetscMemzero(K, sizeof(K));
274 PISM_CHK(ierr, "PetscMemzero");
275
276 for (unsigned int q=0; q<Nq; q++) {
277 for (unsigned int k = 0; k < Nk; k++) { // Test functions
278 const fem::Germ &test_qk=test[q][k];
279 for (unsigned int l = 0; l < Nk; l++) { // Trial functions
280 const fem::Germ &test_ql=test[q][l];
281 K[k][l] += W[q]*(m_cL2*test_qk.val*test_ql.val
282 + m_cH1*(test_qk.dx*test_ql.dx + test_qk.dy*test_ql.dy));
283 } // l
284 } // k
285 } // q
286 m_element.add_contribution(&K[0][0], form);
287 } // j
288 } // i
289
290 if (zeroLocs) {
291 zeroLocs.fix_jacobian(form);
292 }
293
294 ierr = MatAssemblyBegin(form, MAT_FINAL_ASSEMBLY);
295 PISM_CHK(ierr, "MatAssemblyBegin");
296
297 ierr = MatAssemblyEnd(form, MAT_FINAL_ASSEMBLY);
298 PISM_CHK(ierr, "MatAssemblyEnd");
299 }
300
301 } // end of namespace inverse
302 } // end of namespace pism