Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
IPMeanSquareFunctional.cc
Go to the documentation of this file.
1// Copyright (C) 2012, 2014, 2015, 2016, 2017, 2020, 2022, 2023 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 "pism/inverse/functional/IPMeanSquareFunctional.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
25namespace pism {
26namespace inverse {
27
28//! Implicitly set the normalization constant for the functional.
29/*! The normalization constant is selected so that if an input
30array::Vector has component vectors all of length \a scale, then the funtional value will be 1. I.e.
31\f[
32c_N^{-1} = \sum_{i} w_i {\tt scale}^2.
33\f]*/
35
36 // The local value of the weights
37 double value = 0;
38
39 if (m_weights) {
41
42 for (auto p = m_grid->points(); p; p.next()) {
43 const int i = p.i(), j = p.j();
44
45 value += (*m_weights)(i, j);
46 }
47 } else {
48 for (auto p = m_grid->points(); p; p.next()) {
49 value += 1;
50 }
51 }
52
53 m_normalization = GlobalSum(m_grid->com, value);
54 m_normalization *= (scale*scale);
55}
56
58
59 // The value of the objective
60 double value = 0;
61
62 array::AccessScope list{&x};
63
64 if (m_weights) {
65 list.add(*m_weights);
66 for (auto p = m_grid->points(); p; p.next()) {
67 const int i = p.i(), j = p.j();
68
69 Vector2d &x_ij = x(i, j);
70 value += (x_ij.u*x_ij.u + x_ij.v*x_ij.v)*(*m_weights)(i, j);
71 }
72 } else {
73 for (auto p = m_grid->points(); p; p.next()) {
74 const int i = p.i(), j = p.j();
75
76 Vector2d &x_ij = x(i, j);
77 value += (x_ij.u*x_ij.u + x_ij.v*x_ij.v);
78 }
79 }
80 value /= m_normalization;
81
82 GlobalSum( m_grid->com, &value, OUTPUT, 1);
83}
84
86
87 // The value of the objective
88 double value = 0;
89
90 array::AccessScope list{&a, &b};
91
92 if (m_weights) {
93 list.add(*m_weights);
94 for (auto p = m_grid->points(); p; p.next()) {
95 const int i = p.i(), j = p.j();
96
97 Vector2d &a_ij = a(i, j);
98 Vector2d &b_ij = b(i, j);
99 value += (a_ij.u*b_ij.u + a_ij.v*b_ij.v)*(*m_weights)(i, j);
100 }
101 } else {
102 for (auto p = m_grid->points(); p; p.next()) {
103 const int i = p.i(), j = p.j();
104
105 Vector2d &a_ij = a(i, j);
106 Vector2d &b_ij = b(i, j);
107 value += (a_ij.u*b_ij.u + a_ij.v*b_ij.v);
108 }
109 }
110 value /= m_normalization;
111
112 GlobalSum( m_grid->com, &value, OUTPUT, 1);
113}
114
116 gradient.set(0);
117
118 array::AccessScope list{&x, &gradient};
119
120 if (m_weights) {
121 list.add(*m_weights);
122 for (auto p = m_grid->points(); p; p.next()) {
123 const int i = p.i(), j = p.j();
124
125 gradient(i, j).u = 2*x(i, j).u*(*m_weights)(i, j) / m_normalization;
126 gradient(i, j).v = 2*x(i, j).v*(*m_weights)(i, j) / m_normalization;
127 }
128 } else {
129 for (auto p = m_grid->points(); p; p.next()) {
130 const int i = p.i(), j = p.j();
131
132 gradient(i, j).u = 2*x(i, j).u / m_normalization;
133 gradient(i, j).v = 2*x(i, j).v / m_normalization;
134 }
135 }
136}
137
138//! Implicitly set the normalization constant for the functional.
139/*! The normalization constant is selected so that if an input
140array::Scalar has entries all equal to \a scale, then the funtional value will be 1. I.e.
141\f[
142c_N^{-1} = \sum_{i} w_i {\tt scale}^2.
143\f]*/
145
146 // The local value of the weights
147 double value = 0;
148
149 if (m_weights) {
151 for (auto p = m_grid->points(); p; p.next()) {
152 const int i = p.i(), j = p.j();
153
154 value += (*m_weights)(i, j);
155 }
156 } else {
157 for (auto p = m_grid->points(); p; p.next()) {
158 value += 1;
159 }
160 }
161
162 m_normalization = GlobalSum(m_grid->com, value);
163 m_normalization *= (scale*scale);
164}
165
167
168 // The value of the objective
169 double value = 0;
170
171 array::AccessScope list(x);
172
173 if (m_weights) {
174 list.add(*m_weights);
175 for (auto p = m_grid->points(); p; p.next()) {
176 const int i = p.i(), j = p.j();
177
178 double &x_ij = x(i, j);
179 value += x_ij*x_ij*(*m_weights)(i, j);
180 }
181 } else {
182 for (auto p = m_grid->points(); p; p.next()) {
183 const int i = p.i(), j = p.j();
184
185 double &x_ij = x(i, j);
186 value += x_ij*x_ij;
187 }
188 }
189 value /= m_normalization;
190
191 GlobalSum(m_grid->com, &value, OUTPUT, 1);
192}
193
195
196 // The value of the objective
197 double value = 0;
198
199 array::AccessScope list{&a, &b};
200
201 if (m_weights) {
202 list.add(*m_weights);
203 for (auto p = m_grid->points(); p; p.next()) {
204 const int i = p.i(), j = p.j();
205
206 value += (a(i, j)*b(i, j))*(*m_weights)(i, j);
207 }
208 } else {
209 for (auto p = m_grid->points(); p; p.next()) {
210 const int i = p.i(), j = p.j();
211
212 value += (a(i, j)*b(i, j));
213 }
214 }
215 value /= m_normalization;
216
217 GlobalSum(m_grid->com, &value, OUTPUT, 1);
218}
219
220
222 gradient.set(0);
223
224 array::AccessScope list{&x, &gradient};
225
226 if (m_weights) {
227 list.add(*m_weights);
228 for (auto p = m_grid->points(); p; p.next()) {
229 const int i = p.i(), j = p.j();
230
231 gradient(i, j) = 2*x(i, j)*(*m_weights)(i, j) / m_normalization;
232 }
233 } else {
234 for (auto p = m_grid->points(); p; p.next()) {
235 const int i = p.i(), j = p.j();
236
237 gradient(i, j) = 2*x(i, j) / m_normalization;
238 }
239 }
240}
241
242} // end of namespace inverse
243} // end of namespace pism
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition Vector2d.hh:29
void add(const PetscAccessible &v)
Definition Array.cc:896
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void add(double alpha, const Array2D< T > &x)
Definition Array2D.hh:65
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
std::shared_ptr< const Grid > m_grid
virtual void normalize(double scale)
Implicitly set the normalization constant for the functional.
virtual void dot(array::Scalar &a, array::Scalar &b, double *OUTPUT)
Computes the inner product .
virtual void valueAt(array::Scalar &x, double *OUTPUT)
virtual void gradientAt(array::Scalar &x, array::Scalar &gradient)
virtual void gradientAt(array::Vector &x, array::Vector &gradient)
virtual void dot(array::Vector &a, array::Vector &b, double *OUTPUT)
Computes the inner product .
virtual void valueAt(array::Vector &x, double *OUTPUT)
virtual void normalize(double scale)
Implicitly set the normalization constant for the functional.
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)