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