PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
IP_L2NormFunctional.cc
Go to the documentation of this file.
1 // Copyright (C) 2012, 2014, 2015, 2016, 2017, 2020, 2022, 2023 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 "pism/inverse/functional/IP_L2NormFunctional.hh"
20 #include "pism/util/Grid.hh"
21 #include "pism/util/array/Vector.hh"
22 #include "pism/util/array/Scalar.hh"
23 #include "pism/util/pism_utilities.hh"
24 
25 namespace pism {
26 namespace inverse {
27 
29 
30  const unsigned int Nq = m_element.n_pts();
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_q[Nq_max];
37 
38  array::AccessScope list(x);
39 
40  // Loop through all LOCAL elements.
41  const int
42  xs = m_element_index.lxs,
43  xm = m_element_index.lxm,
44  ys = m_element_index.lys,
45  ym = m_element_index.lym;
46 
47  for (int j = ys; j < ys + ym; j++) {
48  for (int i = xs; i < xs + xm; i++) {
49  m_element.reset(i, j);
50 
51  // Obtain values of x at the quadrature points for the element.
52  double tmp[fem::q1::n_chi];
53  m_element.nodal_values(x.array(), tmp);
54  m_element.evaluate(tmp, x_q);
55 
56  for (unsigned int q = 0; q < Nq; q++) {
57  auto W = m_element.weight(q);
58  const double x_qq = x_q[q];
59  value += W*x_qq*x_qq;
60  } // q
61  } // j
62  } // i
63 
64  GlobalSum(m_grid->com, &value, OUTPUT, 1);
65 }
66 
68 
69  const unsigned int Nq = m_element.n_pts();
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  array::AccessScope list{&a, &b};
79 
80  // Loop through all LOCAL elements.
81  const int
82  xs = m_element_index.lxs,
83  xm = m_element_index.lxm,
84  ys = m_element_index.lys,
85  ym = m_element_index.lym;
86 
87  for (int j = ys; j < ys + ym; j++) {
88  for (int i = xs; i < xs + xm; i++) {
89  m_element.reset(i, j);
90 
91  double tmp[fem::q1::n_chi];
92  m_element.nodal_values(a.array(), tmp);
93  m_element.evaluate(tmp, a_q);
94 
95  m_element.nodal_values(b.array(), tmp);
96  m_element.evaluate(tmp, b_q);
97 
98  for (unsigned int q = 0; q < Nq; q++) {
99  auto W = m_element.weight(q);
100  value += W*a_q[q]*b_q[q];
101  } // q
102  } // j
103  } // i
104 
105  GlobalSum(m_grid->com, &value, OUTPUT, 1);
106 }
107 
109 
110  const unsigned int Nk = fem::q1::n_chi;
111  const unsigned int Nq = m_element.n_pts();
112  const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
113 
114  // Clear the gradient before doing anything with it!
115  gradient.set(0);
116 
117  double x_q[Nq_max];
118  double gradient_e[Nk];
119 
120  array::AccessScope list{&x, &gradient};
121 
122  // Loop through all local and ghosted elements.
123  const int
124  xs = m_element_index.xs,
125  xm = m_element_index.xm,
126  ys = m_element_index.ys,
127  ym = m_element_index.ym;
128 
129  for (int j = ys; j < ys + ym; j++) {
130  for (int i = xs; i < xs + xm; i++) {
131 
132  // Reset the DOF map for this element.
133  m_element.reset(i, j);
134 
135  // Obtain values of x at the quadrature points for the element.
136  double tmp[Nk];
137  m_element.nodal_values(x.array(), tmp);
138  m_element.evaluate(tmp, x_q);
139 
140  // Zero out the element-local residual in prep for updating it.
141  for (unsigned int k = 0; k < Nk; k++) {
142  gradient_e[k] = 0;
143  }
144 
145  for (unsigned int q = 0; q < Nq; q++) {
146  auto W = m_element.weight(q);
147  const double x_qq = x_q[q];
148  for (unsigned int k = 0; k < Nk; k++) {
149  gradient_e[k] += 2*W*x_qq*m_element.chi(q, k).val;
150  } // k
151  } // q
152  m_element.add_contribution(gradient_e, gradient.array());
153  } // j
154  } // i
155 }
156 
158 
159  const unsigned int Nq = m_element.n_pts();
160  const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
161 
162  // The value of the objective
163  double value = 0;
164 
165  Vector2d x_q[Nq_max];
166 
167  array::AccessScope list(x);
168 
169  // Loop through all local and ghosted elements.
170  const int
171  xs = m_element_index.lxs,
172  xm = m_element_index.lxm,
173  ys = m_element_index.lys,
174  ym = m_element_index.lym;
175 
176  for (int j = ys; j < ys + ym; j++) {
177  for (int i = xs; i < xs + xm; i++) {
178  m_element.reset(i, j);
179 
180  // Obtain values of x at the quadrature points for the element.
182  m_element.nodal_values(x.array(), tmp);
183  m_element.evaluate(tmp, x_q);
184 
185  for (unsigned int q = 0; q < Nq; q++) {
186  auto W = m_element.weight(q);
187  const Vector2d &x_qq = x_q[q];
188  value += W*(x_qq.u*x_qq.u + x_qq.v*x_qq.v);
189  } // q
190  } // j
191  } // i
192 
193  GlobalSum(m_grid->com, &value, OUTPUT, 1);
194 }
195 
197 
198  const unsigned int Nq = m_element.n_pts();
199  const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
200 
201  // The value of the objective
202  double value = 0;
203 
204  Vector2d a_q[Nq_max];
205  Vector2d b_q[Nq_max];
206 
207  array::AccessScope list{&a, &b};
208 
209  // Loop through all LOCAL elements.
210  const int
211  xs = m_element_index.lxs,
212  xm = m_element_index.lxm,
213  ys = m_element_index.lys,
214  ym = m_element_index.lym;
215 
216  for (int j = ys; j < ys + ym; j++) {
217  for (int i = xs; i < xs + xm; i++) {
218  m_element.reset(i, j);
219 
220  // Obtain values of x at the quadrature points for the element.
222  m_element.nodal_values(a.array(), tmp);
223  m_element.evaluate(tmp, a_q);
224  m_element.nodal_values(b.array(), tmp);
225  m_element.evaluate(tmp, b_q);
226 
227  for (unsigned int q = 0; q < Nq; q++) {
228  auto W = m_element.weight(q);
229  value += W*(a_q[q].u*b_q[q].u + a_q[q].v*b_q[q].v);
230  } // q
231  } // j
232  } // i
233 
234  GlobalSum(m_grid->com, &value, OUTPUT, 1);
235 }
236 
238 
239  const unsigned int Nk = fem::q1::n_chi;
240  const unsigned int Nq = m_element.n_pts();
241  const unsigned int Nq_max = fem::MAX_QUADRATURE_SIZE;
242 
243  // Clear the gradient before doing anything with it!
244  gradient.set(0);
245 
246  Vector2d x_q[Nq_max];
247  Vector2d gradient_e[Nk];
248 
249  array::AccessScope list{&x, &gradient};
250 
251  // Loop through all local and ghosted elements.
252  const int
253  xs = m_element_index.xs,
254  xm = m_element_index.xm,
255  ys = m_element_index.ys,
256  ym = m_element_index.ym;
257 
258  for (int j = ys; j < ys + ym; j++) {
259  for (int i = xs; i < xs + xm; i++) {
260 
261  // Reset the DOF map for this element.
262  m_element.reset(i, j);
263 
264  // Obtain values of x at the quadrature points for the element.
265  Vector2d tmp[Nk];
266  m_element.nodal_values(x.array(), tmp);
267  m_element.evaluate(tmp, x_q);
268 
269  // Zero out the element-local residual in prep for updating it.
270  for (unsigned int k = 0; k < Nk; k++) {
271  gradient_e[k].u = 0;
272  gradient_e[k].v = 0;
273  }
274 
275  for (unsigned int q = 0; q < Nq; q++) {
276  auto W = m_element.weight(q);
277  const Vector2d &x_qq = x_q[q];
278  for (unsigned int k = 0; k < Nk; k++) {
279  double gcommon = 2*W*m_element.chi(q, k).val;
280  gradient_e[k].u += gcommon*x_qq.u;
281  gradient_e[k].v += gcommon*x_qq.v;
282  } // k
283  } // q
284  m_element.add_contribution(gradient_e, gradient.array());
285  } // j
286  } // i
287 }
288 
289 } // end of namespace inverse
290 } // end of namespace pism
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2d.hh:29
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void evaluate(const T *x, T *vals, T *dx, T *dy)
Given nodal values, compute the values and partial derivatives at the quadrature points.
Definition: Element.hh:195
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
Definition: Element.cc:187
void add_contribution(const T *local, T **y_global) const
Add the values of element-local contributions y to the global vector y_global.
Definition: Element.hh:231
void nodal_values(const array::Scalar &x_global, int *result) const
Get nodal values of an integer mask.
Definition: Element.cc:176
int xm
total number of elements to loop over in the x-direction.
int lym
total number local elements in y direction.
int lxm
total number local elements in x direction.
int lxs
x-index of the first local element.
int ym
total number of elements to loop over in the y-direction.
int ys
y-coordinate of the first element to loop over.
int lys
y-index of the first local element.
int xs
x-coordinate of the first element to loop over.
const Germ & chi(unsigned int q, unsigned int k) const
Definition: Element.hh:73
int n_pts() const
Number of quadrature points.
Definition: Element.hh:80
double weight(unsigned int q) const
Weight of the quadrature point q
Definition: Element.hh:85
std::shared_ptr< const Grid > m_grid
Definition: IPFunctional.hh:71
virtual void valueAt(array::Scalar &x, double *OUTPUT)
Computes the value of the functional at the vector x.
virtual void gradientAt(array::Scalar &x, array::Scalar &gradient)
Computes the gradient of the functional at the vector x.
virtual void dot(array::Scalar &a, array::Scalar &b, double *v)
Computes the inner product .
virtual void gradientAt(array::Vector &x, array::Vector &gradient)
Computes the gradient of the functional at the vector x.
virtual void valueAt(array::Vector &x, double *v)
Computes the value of the functional at the vector x.
virtual void dot(array::Vector &a, array::Vector &b, double *v)
Computes the inner product .
const int n_chi
Definition: FEM.hh:191
const unsigned int MAX_QUADRATURE_SIZE
Definition: FEM.hh:173
static const double k
Definition: exactTestP.cc:42
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
double val
Function value.
Definition: FEM.hh:153