PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
iceModelVec2.cc
Go to the documentation of this file.
1 // Copyright (C) 2008--2021 Ed Bueler 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 <cstring>
20 #include <cstdlib>
21 
22 #include <cassert>
23 
24 #include <memory>
25 
26 #include "iceModelVec.hh"
27 #include "iceModelVec_helpers.hh"
28 #include "IceModelVec2V.hh"
29 #include "pism/util/IceModelVec_impl.hh"
30 
31 #include "IceGrid.hh"
32 #include "ConfigInterface.hh"
33 
34 #include "error_handling.hh"
35 
36 #include "pism_utilities.hh"
37 #include "io/io_helpers.hh"
38 #include "pism/util/io/File.hh"
39 #include "pism/util/Logger.hh"
40 #include "pism/util/Context.hh"
41 #include "pism/util/VariableMetadata.hh"
42 
43 namespace pism {
44 
45 // this file contains methods for derived classes IceModelVec2S and IceModelVec2Int
46 
47 // methods for base class IceModelVec are in "iceModelVec.cc"
48 
49 IceModelVec2S::IceModelVec2S(IceGrid::ConstPtr grid, const std::string &name,
50  IceModelVecKind ghostedp, int width)
51  : IceModelVec(grid, name, ghostedp, 1, width, {0.0}) {
52  set_begin_access_use_dof(false);
53 }
54 
55 std::shared_ptr<IceModelVec2S> duplicate(const IceModelVec2S &source) {
56  auto result = std::make_shared<IceModelVec2S>(source.grid(),
57  source.get_name(),
59  1);
60  result->metadata() = source.metadata();
61 
62  return result;
63 }
64 
66  return static_cast<double**>(m_array);
67 }
68 
69 double const* const* IceModelVec2S::array() const {
70  return static_cast<double const* const*>(m_array);
71 }
72 
73 
74 //! Sets an IceModelVec2 to the magnitude of a 2D vector field with components `v_x` and `v_y`.
75 /*! Computes the magnitude \b pointwise, so any of v_x, v_y and the IceModelVec
76  this is called on can be the same.
77 
78  Does not communicate.
79  */
81  const IceModelVec2S &v_y,
82  IceModelVec2S &result) {
83  IceModelVec::AccessList list{&result, &v_x, &v_y};
84 
85  for (Points p(*result.grid()); p; p.next()) {
86  const int i = p.i(), j = p.j();
87 
88  result(i, j) = Vector2(v_x(i, j), v_y(i, j)).magnitude();
89  }
90 
91  result.inc_state_counter(); // mark as modified
92 }
93 
94 void compute_magnitude(const IceModelVec2V &input, IceModelVec2S &result) {
95  IceModelVec::AccessList list{&result, &input};
96 
97  for (Points p(*result.grid()); p; p.next()) {
98  const int i = p.i(), j = p.j();
99 
100  result(i, j) = input(i, j).magnitude();
101  }
102 
103  result.inc_state_counter(); // mark as modified
104 }
105 
106 //! Masks out all the areas where \f$ M \le 0 \f$ by setting them to `fill`.
107 void apply_mask(const IceModelVec2S &M, double fill, IceModelVec2S &result) {
108  IceModelVec::AccessList list{&result, &M};
109 
110  for (Points p(*result.grid()); p; p.next()) {
111  const int i = p.i(), j = p.j();
112 
113  if (M(i, j) <= 0.0) {
114  result(i,j) = fill;
115  }
116  }
117 
118  result.inc_state_counter(); // mark as modified
119 }
120 
121 //! \brief Returns the x-derivative at i,j approximated using centered finite
122 //! differences.
123 double diff_x(const IceModelVec2S &array, int i, int j) {
124  return (array(i + 1,j) - array(i - 1,j)) / (2 * array.grid()->dx());
125 }
126 
127 //! \brief Returns the y-derivative at i,j approximated using centered finite
128 //! differences.
129 double diff_y(const IceModelVec2S &array, int i, int j) {
130  return (array(i,j + 1) - array(i,j - 1)) / (2 * array.grid()->dy());
131 }
132 
133 //! \brief Returns the x-derivative at i,j approximated using centered finite
134 //! differences. Respects grid periodicity and uses one-sided FD at grid edges
135 //! if necessary.
136 double diff_x_p(const IceModelVec2S &array, int i, int j) {
137  const auto &grid = *array.grid();
138 
139  if ((grid.periodicity() & X_PERIODIC) != 0) {
140  return diff_x(array, i,j);
141  }
142 
143  if (i == 0) {
144  return (array(i + 1,j) - array(i,j)) / (grid.dx());
145  } else if (i == (int)grid.Mx() - 1) {
146  return (array(i,j) - array(i - 1,j)) / (grid.dx());
147  } else {
148  return diff_x(array, i,j);
149  }
150 }
151 
152 //! \brief Returns the y-derivative at i,j approximated using centered finite
153 //! differences. Respects grid periodicity and uses one-sided FD at grid edges
154 //! if necessary.
155 double diff_y_p(const IceModelVec2S &array, int i, int j) {
156  const auto &grid = *array.grid();
157 
158  if ((grid.periodicity() & Y_PERIODIC) != 0) {
159  return diff_y(array, i,j);
160  }
161 
162  if (j == 0) {
163  return (array(i,j + 1) - array(i,j)) / (grid.dy());
164  }
165 
166  if (j == (int)grid.My() - 1) {
167  return (array(i,j) - array(i,j - 1)) / (grid.dy());
168  }
169 
170  return diff_y(array, i,j);
171 }
172 
173 //! Sums up all the values in an IceModelVec2S object. Ignores ghosts.
174 /*! Avoids copying to a "global" vector.
175  */
176 double sum(const IceModelVec2S &input) {
177  double result = 0;
178 
179  IceModelVec::AccessList list(input);
180 
181  // sum up the local part:
182  for (Points p(*input.grid()); p; p.next()) {
183  result += input(p.i(), p.j());
184  }
185 
186  // find the global sum:
187  return GlobalSum(input.grid()->com, result);
188 }
189 
190 //! Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
191 double max(const IceModelVec2S &input) {
192  IceModelVec::AccessList list(input);
193 
194  auto grid = input.grid();
195 
196  double result = input(grid->xs(), grid->ys());
197  for (Points p(*grid); p; p.next()) {
198  result = std::max(result, input(p.i(), p.j()));
199  }
200 
201  return GlobalMax(grid->com, result);
202 }
203 
204 //! Finds maximum over all the absolute values in an IceModelVec2S object. Ignores ghosts.
205 double absmax(const IceModelVec2S &input) {
206 
207  double result = 0.0;
208 
209  IceModelVec::AccessList list(input);
210  for (Points p(*input.grid()); p; p.next()) {
211  result = std::max(result, std::abs(input(p.i(), p.j())));
212  }
213 
214  return GlobalMax(input.grid()->com, result);
215 }
216 
217 
218 //! Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
219 double min(const IceModelVec2S &input) {
220  IceModelVec::AccessList list(input);
221 
222  auto grid = input.grid();
223 
224  double result = input(grid->xs(), grid->ys());
225  for (Points p(*grid); p; p.next()) {
226  result = std::min(result, input(p.i(), p.j()));
227  }
228 
229  return GlobalMin(grid->com, result);
230 }
231 
232 void IceModelVec2S::add(double alpha, const IceModelVec2S &x) {
233  vec::add(*this, alpha, x, *this);
234 }
235 
236 void IceModelVec2S::add(double alpha, const IceModelVec2S &x, IceModelVec2S &result) const {
237  vec::add(*this, alpha, x, result);
238 }
239 
241  vec::copy(source, *this);
242 }
243 
244 // IceModelVec2Stag
245 
247  IceModelVecKind ghostedp,
248  unsigned int stencil_width)
249  : IceModelVec3(grid, name, ghostedp, 2, stencil_width) {
251 }
252 
253 std::array<double,2> absmax(const IceModelVec2Stag &input) {
254 
255  double z[2] = {0.0, 0.0};
256 
257  IceModelVec::AccessList list(input);
258  for (Points p(*input.grid()); p; p.next()) {
259  const int i = p.i(), j = p.j();
260 
261  z[0] = std::max(z[0], std::abs(input(i, j, 0)));
262  z[1] = std::max(z[1], std::abs(input(i, j, 1)));
263  }
264 
265  double result[2];
266  GlobalMax(input.grid()->com, z, result, 2);
267 
268  return {result[0], result[1]};
269 }
270 
272  IceModelVecKind ghostedp, int width)
273  : IceModelVec2S(grid, name, ghostedp, width) {
275 }
276 
277 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
IceModelVec2Int(IceGrid::ConstPtr grid, const std::string &name, IceModelVecKind ghostedp, int width=1)
void add(double alpha, const IceModelVec2S &x)
IceModelVec2S(IceGrid::ConstPtr grid, const std::string &name, IceModelVecKind ghostedp, int width=1)
Definition: iceModelVec2.cc:49
void copy_from(const IceModelVec2S &source)
IceModelVec2Stag(IceGrid::ConstPtr grid, const std::string &name, IceModelVecKind ghostedp, unsigned int stencil_width=1)
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: iceModelVec.hh:449
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void set_begin_access_use_dof(bool flag)
Definition: iceModelVec.cc:181
const std::string & get_name() const
Get the name of an IceModelVec object.
Definition: iceModelVec.cc:386
void inc_state_counter()
Increment the object state counter.
Definition: iceModelVec.cc:147
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Definition: iceModelVec.hh:202
double magnitude() const
Magnitude.
Definition: Vector2.hh:40
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
void add(const V &x, double alpha, const V &y, V &result, bool scatter=true)
Computes result = x + alpha * y, where x, y, and z are 2D IceModelVecs (scalar or vector).
void copy(const V &source, V &destination, bool scatter=true)
void apply_mask(const IceModelVec2S &M, double fill, IceModelVec2S &result)
Masks out all the areas where by setting them to fill.
double diff_x(const IceModelVec2S &array, int i, int j)
Returns the x-derivative at i,j approximated using centered finite differences.
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
@ Y_PERIODIC
Definition: IceGrid.hh:48
@ X_PERIODIC
Definition: IceGrid.hh:48
void compute_magnitude(const IceModelVec2S &v_x, const IceModelVec2S &v_y, IceModelVec2S &result)
Sets an IceModelVec2 to the magnitude of a 2D vector field with components v_x and v_y.
Definition: iceModelVec2.cc:80
double diff_x_p(const IceModelVec2S &array, int i, int j)
Returns the x-derivative at i,j approximated using centered finite differences. Respects grid periodi...
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
double absmax(const IceModelVec2S &input)
Finds maximum over all the absolute values in an IceModelVec2S object. Ignores ghosts.
std::shared_ptr< IceModelVec2S > duplicate(const IceModelVec2S &source)
Definition: iceModelVec2.cc:55
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
void GlobalMin(MPI_Comm comm, double *local, double *result, int count)
double diff_y_p(const IceModelVec2S &array, int i, int j)
Returns the y-derivative at i,j approximated using centered finite differences. Respects grid periodi...
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
double diff_y(const IceModelVec2S &array, int i, int j)
Returns the y-derivative at i,j approximated using centered finite differences.
IceModelVecKind
What "kind" of a vector to create: with or without ghosts.
Definition: iceModelVec.hh:49
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
double sum(const IceModelVec2S &input)
Sums up all the values in an IceModelVec2S object. Ignores ghosts.
InterpolationType interpolation_type