PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
IceGrid.hh
Go to the documentation of this file.
1 // Copyright (C) 2004-2021 Jed Brown, 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 #ifndef PISM_GRID_H
20 #define PISM_GRID_H
21 
22 #include <cassert>
23 #include <vector>
24 #include <string>
25 #include <memory> // shared_ptr
26 
27 #include <mpi.h> // MPI_Comm
28 
29 namespace pism {
30 
31 class Context;
32 class Config;
33 
34 namespace petsc {
35 class DM;
36 } // end of namespace petsc
37 
38 class File;
39 namespace units {
40 class System;
41 }
42 class Vars;
43 class Logger;
44 
45 class MappingInfo;
46 
47 typedef enum {UNKNOWN = 0, EQUAL, QUADRATIC} SpacingType;
48 typedef enum {NOT_PERIODIC = 0, X_PERIODIC = 1, Y_PERIODIC = 2, XY_PERIODIC = 3} Periodicity;
49 
51 
52 GridRegistration string_to_registration(const std::string &keyword);
53 std::string registration_to_string(GridRegistration registration);
54 
55 Periodicity string_to_periodicity(const std::string &keyword);
56 std::string periodicity_to_string(Periodicity p);
57 
58 SpacingType string_to_spacing(const std::string &keyword);
59 std::string spacing_to_string(SpacingType s);
60 
61 //! @brief Contains parameters of an input file grid.
62 class grid_info {
63 public:
64  grid_info();
65  grid_info(const File &file, const std::string &variable,
66  std::shared_ptr<units::System> unit_system,
67  GridRegistration registration);
68  void report(const Logger &log, int threshold, std::shared_ptr<units::System> s) const;
69  // dimension lengths
70  unsigned int t_len, x_len, y_len, z_len;
71  //! current time (seconds)
72  double time;
73  //! x-coordinate of the domain center
74  double x0;
75  //! y-coordinate of the domain center
76  double y0;
77  //! domain half-width
78  double Lx;
79  //! domain half-height
80  double Ly;
81  //! minimal value of the z dimension
82  double z_min;
83  //! maximal value of the z dimension
84  double z_max;
85 
86  //! x coordinates
87  std::vector<double> x;
88  //! y coordinates
89  std::vector<double> y;
90  //! z coordinates
91  std::vector<double> z;
92  std::string filename;
93 private:
94  void reset();
95 };
96 
97 //! Grid parameters; used to collect defaults before an IceGrid is allocated.
98 /* Make sure that all of
99  - `horizontal_size_from_options()`
100  - `horizontal_extent_from_options()`
101  - `vertical_grid_from_options()`
102  - `ownership_ranges_from_options()`
103 
104  are called *or* all data members (`Lx`, `Ly`, `x0`, `y0`, `Mx`, `My`, `z`, `periodicity`,
105  `procs_x`, `procs_y`) are set manually before using an instance of GridParameters.
106 
107  Call `validate()` to check current parameters.
108 */
110 public:
111  //! Create an uninitialized GridParameters instance.
112  GridParameters();
113 
114  //! Initialize grid defaults from a configuration database.
115  GridParameters(std::shared_ptr<const Config> config);
116 
117  //! Initialize grid defaults from a configuration database and a NetCDF variable.
118  GridParameters(std::shared_ptr<const Context> ctx,
119  const std::string &filename,
120  const std::string &variable_name,
121  GridRegistration r);
122  //! Initialize grid defaults from a configuration database and a NetCDF variable.
123  GridParameters(std::shared_ptr<const Context> ctx,
124  const File &file,
125  const std::string &variable_name,
126  GridRegistration r);
127 
128  //! Process -Mx and -My; set Mx and My.
130  //! Process -Lx, -Ly, -x0, -y0, -x_range, -y_range; set Lx, Ly, x0, y0.
131  void horizontal_extent_from_options(std::shared_ptr<units::System> unit_system);
132  //! Process -Mz and -Lz; set z;
133  void vertical_grid_from_options(std::shared_ptr<const Config> config);
134  //! Re-compute ownership ranges. Uses current values of Mx and My.
135  void ownership_ranges_from_options(unsigned int size);
136 
137  //! Validate data members.
138  void validate() const;
139 
140  //! Domain half-width in the X direction.
141  double Lx;
142  //! Domain half-width in the Y direction.
143  double Ly;
144  //! Domain center in the X direction.
145  double x0;
146  //! Domain center in the Y direction.
147  double y0;
148  //! Number of grid points in the X direction.
149  unsigned int Mx;
150  //! Number of grid points in the Y direction.
151  unsigned int My;
152  //! Grid registration.
154  //! Grid periodicity.
156  //! Vertical levels.
157  std::vector<double> z;
158  //! Processor ownership ranges in the X direction.
159  std::vector<unsigned int> procs_x;
160  //! Processor ownership ranges in the Y direction.
161  std::vector<unsigned int> procs_y;
162 private:
163  void init_from_config(std::shared_ptr<const Config> config);
164  void init_from_file(std::shared_ptr<const Context> ctx, const File &file,
165  const std::string &variable_name,
166  GridRegistration r);
167 };
168 
169 //! Describes the PISM grid and the distribution of data across processors.
170 /*!
171  This class holds parameters describing the grid, including the vertical
172  spacing and which part of the horizontal grid is owned by the processor. It
173  contains the dimensions of the PISM (4-dimensional, x*y*z*time) computational
174  box. The vertical spacing can be quite arbitrary.
175 
176  It creates and destroys a two dimensional `PETSc` `DA` (distributed array).
177  The creation of this `DA` is the point at which PISM gets distributed across
178  multiple processors.
179 
180  \section computational_grid Organization of PISM's computational grid
181 
182  PISM uses the class IceGrid to manage computational grids and their
183  parameters.
184 
185  Computational grids PISM can use are
186  - rectangular,
187  - equally spaced in the horizintal (X and Y) directions,
188  - distributed across processors in horizontal dimensions only (every column
189  is stored on one processor only),
190  - are periodic in both X and Y directions (in the topological sence).
191 
192  Each processor "owns" a rectangular patch of `xm` times `ym` grid points with
193  indices starting at `xs` and `ys` in the X and Y directions respectively.
194 
195  The typical code performing a point-wise computation will look like
196 
197  \code
198  for (Points p(grid); p; p.next()) {
199  const int i = p.i(), j = p.j();
200  field(i,j) = value;
201  }
202  \endcode
203 
204  For finite difference (and some other) computations we often need to know
205  values at map-plane neighbors of a grid point.
206 
207  We say that a patch owned by a processor is surrounded by a strip of "ghost"
208  grid points belonging to patches next to the one in question. This lets us to
209  access (read) values at all the eight neighbors of a grid point for *all*
210  the grid points, including ones at an edge of a processor patch *and* at an
211  edge of a computational domain.
212 
213  All the values *written* to ghost points will be lost next time ghost values
214  are updated.
215 
216  Sometimes it is beneficial to update ghost values locally (for instance when
217  a computation A uses finite differences to compute derivatives of a quantity
218  produced using a purely local (point-wise) computation B). In this case the
219  loop above can be modified to look like
220 
221  \code
222  for (PointsWithGhosts p(grid, ghost_width); p; p.next()) {
223  const int i = p.i(), j = p.j();
224  field(i,j) = value;
225  }
226  \endcode
227 */
228 class IceGrid {
229 public:
230  ~IceGrid();
231 
232  typedef std::shared_ptr<IceGrid> Ptr;
233  typedef std::shared_ptr<const IceGrid> ConstPtr;
234 
235  IceGrid(std::shared_ptr<const Context> ctx, const GridParameters &p);
236 
237  static std::vector<double> compute_vertical_levels(double new_Lz, unsigned int new_Mz,
238  SpacingType spacing, double Lambda = 0.0);
239 
240  static Ptr Shallow(std::shared_ptr<const Context> ctx,
241  double Lx, double Ly,
242  double x0, double y0,
243  unsigned int Mx, unsigned int My,
245  Periodicity p);
246 
247  static Ptr FromFile(std::shared_ptr<const Context> ctx,
248  const File &file, const std::string &var_name,
249  GridRegistration r);
250 
251  static Ptr FromFile(std::shared_ptr<const Context> ctx,
252  const std::string &file, const std::vector<std::string> &var_names,
253  GridRegistration r);
254 
255  static Ptr FromOptions(std::shared_ptr<const Context> ctx);
256 
257  std::shared_ptr<petsc::DM> get_dm(int dm_dof, int stencil_width) const;
258 
259  void report_parameters() const;
260 
261  void compute_point_neighbors(double X, double Y,
262  int &i_left, int &i_right,
263  int &j_bottom, int &j_top) const;
264  std::vector<int> compute_point_neighbors(double X, double Y) const;
265  std::vector<double> compute_interp_weights(double x, double y) const;
266 
267  unsigned int kBelowHeight(double height) const;
268 
269  std::shared_ptr<const Context> ctx() const;
270 
271  int xs() const;
272  int xm() const;
273  int ys() const;
274  int ym() const;
275 
276  const std::vector<double>& x() const;
277  double x(size_t i) const;
278 
279  const std::vector<double>& y() const;
280  double y(size_t i) const;
281 
282  const std::vector<double>& z() const;
283  double z(size_t i) const;
284 
285  double dx() const;
286  double dy() const;
287  double cell_area() const;
288 
289  unsigned int Mx() const;
290  unsigned int My() const;
291  unsigned int Mz() const;
292 
293  double Lx() const;
294  double Ly() const;
295  double Lz() const;
296  double x0() const;
297  double y0() const;
298 
299  const MappingInfo& get_mapping_info() const;
300  void set_mapping_info(const MappingInfo &info);
301 
302  double dz_min() const;
303  double dz_max() const;
304 
305  Periodicity periodicity() const;
307 
308  unsigned int size() const;
309  int rank() const;
310 
311  const MPI_Comm com;
312 
313  Vars& variables();
314  const Vars& variables() const;
315 
316  int pio_io_decomposition(int dof, int output_datatype) const;
317 
318  //! Maximum number of degrees of freedom supported by PISM.
319  /*!
320  * This is also the maximum number of records an IceModelVec2T can hold.
321  */
322  static const int max_dm_dof = 10000;
323 
324  //! Maximum stencil width supported by PISM.
325  static const int max_stencil_width = 10000;
326 private:
327  struct Impl;
329 
330  // Hide copy constructor / assignment operator.
331  IceGrid(const IceGrid &);
333 };
334 
335 double radius(const IceGrid &grid, int i, int j);
336 
337 //! @brief Check if a point `(i,j)` is in the strip of `stripwidth`
338 //! meters around the edge of the computational domain.
339 inline bool in_null_strip(const IceGrid& grid, int i, int j, double strip_width) {
340  return (strip_width > 0.0 &&
341  (grid.x(i) <= grid.x(0) + strip_width ||
342  grid.x(i) >= grid.x(grid.Mx()-1) - strip_width ||
343  grid.y(j) <= grid.y(0) + strip_width ||
344  grid.y(j) >= grid.y(grid.My()-1) - strip_width));
345 }
346 
347 /*!
348  * Return `true` if a point `(i, j)` is at an edge of `grid`.
349  */
350 inline bool grid_edge(const IceGrid &grid, int i, int j) {
351  return ((j == 0) or (j == (int)grid.My() - 1) or
352  (i == 0) or (i == (int)grid.Mx() - 1));
353 }
354 
355 /** Iterator class for traversing the grid, including ghost points.
356  *
357  * Usage:
358  *
359  * `for (PointsWithGhosts p(grid, stencil_width); p; p.next()) { ... }`
360  */
362 public:
363  PointsWithGhosts(const IceGrid &g, unsigned int stencil_width = 1) {
364  m_i_first = g.xs() - stencil_width;
365  m_i_last = g.xs() + g.xm() + stencil_width - 1;
366  m_j_first = g.ys() - stencil_width;
367  m_j_last = g.ys() + g.ym() + stencil_width - 1;
368 
369  m_i = m_i_first;
370  m_j = m_j_first;
371  m_done = false;
372  }
373 
374  int i() const {
375  return m_i;
376  }
377  int j() const {
378  return m_j;
379  }
380 
381  void next() {
382  assert(not m_done);
383  m_i += 1;
384  if (m_i > m_i_last) {
385  m_i = m_i_first; // wrap around
386  m_j += 1;
387  }
388  if (m_j > m_j_last) {
389  m_j = m_j_first; // ensure that indexes are valid
390  m_done = true;
391  }
392  }
393 
394  operator bool() const {
395  return not m_done;
396  }
397 private:
398  int m_i, m_j;
400  bool m_done;
401 };
402 
403 /** Iterator class for traversing the grid (without ghost points).
404  *
405  * Usage:
406  *
407  * `for (Points p(grid); p; p.next()) { double foo = p.i(); ... }`
408  */
409 class Points : public PointsWithGhosts {
410 public:
411  Points(const IceGrid &g) : PointsWithGhosts(g, 0) {}
412 };
413 
414 } // end of namespace pism
415 
416 #endif /* PISM_GRID_H */
High-level PISM I/O class.
Definition: File.hh:51
void init_from_file(std::shared_ptr< const Context > ctx, const File &file, const std::string &variable_name, GridRegistration r)
Definition: IceGrid.cc:1229
double x0
Domain center in the X direction.
Definition: IceGrid.hh:145
void ownership_ranges_from_options(unsigned int size)
Re-compute ownership ranges. Uses current values of Mx and My.
Definition: IceGrid.cc:1201
std::vector< unsigned int > procs_y
Processor ownership ranges in the Y direction.
Definition: IceGrid.hh:161
void horizontal_extent_from_options(std::shared_ptr< units::System > unit_system)
Process -Lx, -Ly, -x0, -y0, -x_range, -y_range; set Lx, Ly, x0, y0.
Definition: IceGrid.cc:1272
double Ly
Domain half-width in the Y direction.
Definition: IceGrid.hh:143
Periodicity periodicity
Grid periodicity.
Definition: IceGrid.hh:155
std::vector< double > z
Vertical levels.
Definition: IceGrid.hh:157
GridRegistration registration
Grid registration.
Definition: IceGrid.hh:153
double Lx
Domain half-width in the X direction.
Definition: IceGrid.hh:141
GridParameters(std::shared_ptr< const Config > config)
Initialize grid defaults from a configuration database.
void init_from_config(std::shared_ptr< const Config > config)
Initialize from a configuration database. Does not try to compute ownership ranges.
Definition: IceGrid.cc:1208
GridParameters()
Create an uninitialized GridParameters instance.
Definition: IceGrid.cc:1181
double y0
Domain center in the Y direction.
Definition: IceGrid.hh:147
void vertical_grid_from_options(std::shared_ptr< const Config > config)
Process -Mz and -Lz; set z;.
Definition: IceGrid.cc:1299
unsigned int Mx
Number of grid points in the X direction.
Definition: IceGrid.hh:149
void validate() const
Validate data members.
Definition: IceGrid.cc:1308
std::vector< unsigned int > procs_x
Processor ownership ranges in the X direction.
Definition: IceGrid.hh:159
void horizontal_size_from_options()
Process -Mx and -My; set Mx and My.
Definition: IceGrid.cc:1267
unsigned int My
Number of grid points in the Y direction.
Definition: IceGrid.hh:151
Grid parameters; used to collect defaults before an IceGrid is allocated.
Definition: IceGrid.hh:109
IceGrid & operator=(const IceGrid &)
const std::vector< double > & y() const
Y-coordinates.
Definition: IceGrid.cc:963
IceGrid(std::shared_ptr< const Context > ctx, const GridParameters &p)
Create a PISM distributed computational grid.
Definition: IceGrid.cc:246
double y0() const
Y-coordinate of the center of the domain.
Definition: IceGrid.cc:1037
GridRegistration registration() const
Definition: IceGrid.cc:861
Vars & variables()
Dictionary of variables (2D and 3D fields) associated with this grid.
Definition: IceGrid.cc:908
int rank() const
MPI rank.
Definition: IceGrid.cc:898
double cell_area() const
Definition: IceGrid.cc:992
Periodicity periodicity() const
Return grid periodicity.
Definition: IceGrid.cc:857
double Lx() const
Half-width of the computational domain.
Definition: IceGrid.cc:1017
static const int max_dm_dof
Maximum number of degrees of freedom supported by PISM.
Definition: IceGrid.hh:322
const std::vector< double > & x() const
X-coordinates.
Definition: IceGrid.cc:953
int ys() const
Global starting index of this processor's subset.
Definition: IceGrid.cc:923
Impl * m_impl
Definition: IceGrid.hh:327
std::shared_ptr< IceGrid > Ptr
Definition: IceGrid.hh:232
unsigned int Mx() const
Total grid size in the X direction.
Definition: IceGrid.cc:938
static Ptr FromOptions(std::shared_ptr< const Context > ctx)
Create a grid using command-line options and (possibly) an input file.
Definition: IceGrid.cc:1349
std::vector< double > compute_interp_weights(double x, double y) const
Compute 4 interpolation weights necessary for linear interpolation from the current grid....
Definition: IceGrid.cc:801
double x0() const
X-coordinate of the center of the domain.
Definition: IceGrid.cc:1032
unsigned int Mz() const
Number of vertical levels.
Definition: IceGrid.cc:948
static Ptr Shallow(std::shared_ptr< const Context > ctx, double Lx, double Ly, double x0, double y0, unsigned int Mx, unsigned int My, GridRegistration r, Periodicity p)
Initialize a uniform, shallow (3 z-levels) grid with half-widths (Lx,Ly) and Mx by My nodes.
Definition: IceGrid.cc:213
int xm() const
Width of this processor's sub-domain.
Definition: IceGrid.cc:928
void set_mapping_info(const MappingInfo &info)
Definition: IceGrid.cc:1435
static const int max_stencil_width
Maximum stencil width supported by PISM.
Definition: IceGrid.hh:325
std::shared_ptr< const Context > ctx() const
Return execution context this grid corresponds to.
Definition: IceGrid.cc:866
double dz_max() const
Maximum vertical spacing.
Definition: IceGrid.cc:1007
IceGrid(const IceGrid &)
int pio_io_decomposition(int dof, int output_datatype) const
Definition: IceGrid.cc:1459
int xs() const
Global starting index of this processor's subset.
Definition: IceGrid.cc:918
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
unsigned int kBelowHeight(double height) const
Return the index k into zlevels[] so that zlevels[k] <= height < zlevels[k+1] and k < Mz.
Definition: IceGrid.cc:435
double dx() const
Horizontal grid spacing.
Definition: IceGrid.cc:983
double dy() const
Horizontal grid spacing.
Definition: IceGrid.cc:988
const std::vector< double > & z() const
Z-coordinates within the ice.
Definition: IceGrid.cc:973
const MPI_Comm com
Definition: IceGrid.hh:311
void compute_point_neighbors(double X, double Y, int &i_left, int &i_right, int &j_bottom, int &j_top) const
Computes indices of grid points to the lower left and upper right from (X,Y).
Definition: IceGrid.cc:770
const MappingInfo & get_mapping_info() const
Definition: IceGrid.cc:1431
static Ptr FromFile(std::shared_ptr< const Context > ctx, const File &file, const std::string &var_name, GridRegistration r)
Create a grid from a file, get information from variable var_name.
Definition: IceGrid.cc:325
static std::vector< double > compute_vertical_levels(double new_Lz, unsigned int new_Mz, SpacingType spacing, double Lambda=0.0)
Set the vertical levels in the ice according to values in Mz (number of levels), Lz (domain height),...
Definition: IceGrid.cc:389
std::shared_ptr< petsc::DM > get_dm(int dm_dof, int stencil_width) const
Get a PETSc DM ("distributed array manager") object for given dof (number of degrees of freedom per g...
Definition: IceGrid.cc:841
void report_parameters() const
Report grid parameters.
Definition: IceGrid.cc:671
double Lz() const
Height of the computational domain.
Definition: IceGrid.cc:1027
double dz_min() const
Minimum vertical spacing.
Definition: IceGrid.cc:997
int ym() const
Width of this processor's sub-domain.
Definition: IceGrid.cc:933
unsigned int size() const
MPI communicator size.
Definition: IceGrid.cc:903
unsigned int My() const
Total grid size in the Y direction.
Definition: IceGrid.cc:943
double Ly() const
Half-width of the computational domain.
Definition: IceGrid.cc:1022
Describes the PISM grid and the distribution of data across processors.
Definition: IceGrid.hh:228
A basic logging class.
Definition: Logger.hh:40
PointsWithGhosts(const IceGrid &g, unsigned int stencil_width=1)
Definition: IceGrid.hh:363
Points(const IceGrid &g)
Definition: IceGrid.hh:411
A class for passing PISM variables from the core to other parts of the code (such as climate couplers...
Definition: Vars.hh:38
double Ly
domain half-height
Definition: IceGrid.hh:80
unsigned int y_len
Definition: IceGrid.hh:70
double Lx
domain half-width
Definition: IceGrid.hh:78
std::string filename
Definition: IceGrid.hh:92
unsigned int x_len
Definition: IceGrid.hh:70
unsigned int z_len
Definition: IceGrid.hh:70
double z_min
minimal value of the z dimension
Definition: IceGrid.hh:82
double time
current time (seconds)
Definition: IceGrid.hh:72
grid_info(const File &file, const std::string &variable, std::shared_ptr< units::System > unit_system, GridRegistration registration)
std::vector< double > z
z coordinates
Definition: IceGrid.hh:91
void report(const Logger &log, int threshold, std::shared_ptr< units::System > s) const
Definition: IceGrid.cc:1072
std::vector< double > x
x coordinates
Definition: IceGrid.hh:87
double x0
x-coordinate of the domain center
Definition: IceGrid.hh:74
double y0
y-coordinate of the domain center
Definition: IceGrid.hh:76
double z_max
maximal value of the z dimension
Definition: IceGrid.hh:84
unsigned int t_len
Definition: IceGrid.hh:70
std::vector< double > y
y coordinates
Definition: IceGrid.hh:89
Contains parameters of an input file grid.
Definition: IceGrid.hh:62
std::string periodicity_to_string(Periodicity p)
Convert Periodicity to a STL string.
Definition: IceGrid.cc:152
double radius(const IceGrid &grid, int i, int j)
Returns the distance from the point (i,j) to the origin.
Definition: IceGrid.cc:1042
Periodicity
Definition: IceGrid.hh:48
@ NOT_PERIODIC
Definition: IceGrid.hh:48
@ XY_PERIODIC
Definition: IceGrid.hh:48
@ Y_PERIODIC
Definition: IceGrid.hh:48
@ X_PERIODIC
Definition: IceGrid.hh:48
static const double g
Definition: exactTestP.cc:39
SpacingType
Definition: IceGrid.hh:47
@ UNKNOWN
Definition: IceGrid.hh:47
@ QUADRATIC
Definition: IceGrid.hh:47
@ EQUAL
Definition: IceGrid.hh:47
std::string spacing_to_string(SpacingType s)
Convert SpacingType to an STL string.
Definition: IceGrid.cc:179
Periodicity string_to_periodicity(const std::string &keyword)
Convert a string to Periodicity.
Definition: IceGrid.cc:129
std::string registration_to_string(GridRegistration registration)
Definition: IceGrid.cc:200
SpacingType string_to_spacing(const std::string &keyword)
Convert an STL string to SpacingType.
Definition: IceGrid.cc:167
GridRegistration
Definition: IceGrid.hh:50
@ CELL_CENTER
Definition: IceGrid.hh:50
@ CELL_CORNER
Definition: IceGrid.hh:50
bool in_null_strip(const IceGrid &grid, int i, int j, double strip_width)
Check if a point (i,j) is in the strip of stripwidth meters around the edge of the computational doma...
Definition: IceGrid.hh:339
GridRegistration string_to_registration(const std::string &keyword)
Definition: IceGrid.cc:189
bool grid_edge(const IceGrid &grid, int i, int j)
Definition: IceGrid.hh:350
Internal structures of IceGrid.
Definition: IceGrid.cc:49