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