Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Grid.hh
Go to the documentation of this file.
1// Copyright (C) 2004-2025 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 "io/File.hh"
23#include <cassert>
24#include <memory> // shared_ptr
25#include <string>
26#include <vector>
27#include <map>
28
29#include <mpi.h> // MPI_Comm
30
31#include "pism/util/Interpolation1D.hh"
32
33namespace pism {
34
35class Config;
36class Context;
37class File;
38class InputInterpolation;
39class Logger;
40class MappingInfo;
41class Vars;
42
43namespace petsc {
44class DM;
45} // end of namespace petsc
46
47namespace units {
48class System;
49}
50
51namespace grid {
52
54typedef enum {NOT_PERIODIC = 0, X_PERIODIC = 1, Y_PERIODIC = 2, XY_PERIODIC = 3} Periodicity;
55
57
58Registration string_to_registration(const std::string &keyword);
59std::string registration_to_string(Registration registration);
60
61Periodicity string_to_periodicity(const std::string &keyword);
63
64VerticalSpacing string_to_spacing(const std::string &keyword);
66
67//! @brief Contains parameters of an input file grid.
69public:
70 InputGridInfo(const File &file, const std::string &variable,
71 std::shared_ptr<units::System> unit_system, Registration registration);
72
73 void report(const Logger &log, int threshold, std::shared_ptr<units::System> s) const;
74
75 // dimension lengths
76 unsigned int t_len;
77 //! x-coordinate of the domain center
78 double x0;
79 //! y-coordinate of the domain center
80 double y0;
81 //! domain half-width
82 double Lx;
83 //! domain half-height
84 double Ly;
85 //! minimal value of the z dimension
86 double z_min;
87 //! maximal value of the z dimension
88 double z_max;
89
90 //! x coordinates
91 std::vector<double> x;
92 //! y coordinates
93 std::vector<double> y;
94 //! z coordinates
95 std::vector<double> z;
96
97 std::string filename;
98
99 /*!
100 * Variable name *found in the file*, which may not match the one expected by PISM.
101 */
102 std::string variable_name;
103
104 std::map<std::string, AxisType> dimension_types;
105
107private:
108 void reset();
109};
110
111//! Grid parameters; used to collect defaults before an Grid is allocated.
112/* Make sure that all of
113 - `horizontal_size_from_options()`
114 - `horizontal_extent_from_options()`
115 - `vertical_grid_from_options()`
116 - `ownership_ranges_from_options()`
117
118 are called *or* all data members (`Lx`, `Ly`, `x0`, `y0`, `Mx`, `My`, `z`, `periodicity`,
119 `procs_x`, `procs_y`) are set manually before using an instance of GridParameters.
120
121 Call `validate()` to check current parameters.
122*/
124public:
125 //! Initialize grid defaults from a configuration database.
126 Parameters(const Config &config);
127
128 //! Initialize grid defaults from a configuration database, but set Mx and My explicitly.
129 Parameters(const Config &config, unsigned Mx_, unsigned My_, double Lx, double Ly);
130
131 //! Initialize grid defaults from a NetCDF variable.
132 Parameters(std::shared_ptr<units::System> unit_system, const File &file,
133 const std::string &variable_name, Registration r);
134
135 static Parameters FromGridDefinition(std::shared_ptr<units::System> unit_system, const File &file,
136 const std::string &variable_name, Registration registration);
137
138 //! Process -Lx, -Ly, -x0, -y0; set Lx, Ly, x0, y0.
140 //! Process -Mz and -Lz; set z;
141 void vertical_grid_from_options(const Config &config);
142 //! Re-compute ownership ranges. Uses current values of Mx and My.
143 void ownership_ranges_from_options(const Config &config, unsigned int size);
144
145 //! Validate data members.
146 void validate() const;
147
148 //! Domain half-width in the X direction
149 double Lx;
150 //! Domain half-width in the Y direction
151 double Ly;
152 //! Domain center in the X direction
153 double x0;
154 //! Domain center in the Y direction
155 double y0;
156 //! Number of grid points in the X direction
157 unsigned int Mx;
158 //! Number of grid points in the Y direction
159 unsigned int My;
160 //! Grid registration
162 //! Grid periodicity
164 //! Vertical levels
165 std::vector<double> z;
166 //! Processor ownership ranges in the X direction
167 std::vector<unsigned int> procs_x;
168 //! Processor ownership ranges in the Y direction
169 std::vector<unsigned int> procs_y;
170
171 //! Name of the variable used to initialize the instance (empty if not used)
172 std::string variable_name;
173private:
174 Parameters() = default;
175};
176} // namespace grid
177
178class Grid;
179
180/** Iterator class for traversing the grid, including ghost points.
181 *
182 * Usage:
183 *
184 * `for (PointsWithGhosts p(grid, stencil_width); p; p.next()) { ... }`
185 */
187public:
188 PointsWithGhosts(const Grid &grid, unsigned int stencil_width = 1);
189
190 int i() const {
191 return m_i;
192 }
193 int j() const {
194 return m_j;
195 }
196
197 void next() {
198 assert(not m_done);
199 m_i += 1;
200 if (m_i > m_i_last) {
201 m_i = m_i_first; // wrap around
202 m_j += 1;
203 }
204 if (m_j > m_j_last) {
205 m_j = m_j_first; // ensure that indexes are valid
206 m_done = true;
207 }
208 }
209
210 operator bool() const {
211 return not m_done;
212 }
213private:
214 int m_i, m_j;
216 bool m_done;
217};
218
219/** Iterator class for traversing the grid (without ghost points).
220 *
221 * Usage:
222 *
223 * `for (Points p(grid); p; p.next()) { int i = p.i(), j = p.j(); ... }`
224 */
225class Points : public PointsWithGhosts {
226public:
227 Points(const Grid &g) : PointsWithGhosts(g, 0) {}
228};
229
230
231//! Describes the PISM grid and the distribution of data across processors.
232/*!
233 This class holds parameters describing the grid, including the vertical
234 spacing and which part of the horizontal grid is owned by the processor. It
235 contains the dimensions of the PISM (4-dimensional, x*y*z*time) computational
236 box. The vertical spacing can be quite arbitrary.
237
238 It creates and destroys a two dimensional `PETSc` `DA` (distributed array).
239 The creation of this `DA` is the point at which PISM gets distributed across
240 multiple processors.
241
242 \section computational_grid Organization of PISM's computational grid
243
244 PISM uses the class Grid to manage computational grids and their
245 parameters.
246
247 Computational grids PISM can use are
248 - rectangular,
249 - equally spaced in the horizintal (X and Y) directions,
250 - distributed across processors in horizontal dimensions only (every column
251 is stored on one processor only),
252 - are periodic in both X and Y directions (in the topological sence).
253
254 Each processor "owns" a rectangular patch of `xm` times `ym` grid points with
255 indices starting at `xs` and `ys` in the X and Y directions respectively.
256
257 The typical code performing a point-wise computation will look like
258
259 \code
260 for (Points p(grid); p; p.next()) {
261 const int i = p.i(), j = p.j();
262 field(i,j) = value;
263 }
264 \endcode
265
266 For finite difference (and some other) computations we often need to know
267 values at map-plane neighbors of a grid point.
268
269 We say that a patch owned by a processor is surrounded by a strip of "ghost"
270 grid points belonging to patches next to the one in question. This lets us to
271 access (read) values at all the eight neighbors of a grid point for *all*
272 the grid points, including ones at an edge of a processor patch *and* at an
273 edge of a computational domain.
274
275 All the values *written* to ghost points will be lost next time ghost values
276 are updated.
277
278 Sometimes it is beneficial to update ghost values locally (for instance when
279 a computation A uses finite differences to compute derivatives of a quantity
280 produced using a purely local (point-wise) computation B). In this case the
281 loop above can be modified to look like
282
283 \code
284 for (PointsWithGhosts p(grid, ghost_width); p; p.next()) {
285 const int i = p.i(), j = p.j();
286 field(i,j) = value;
287 }
288 \endcode
289*/
290class Grid {
291public:
292 ~Grid();
293
294 Grid(std::shared_ptr<const Context> context, const grid::Parameters &p);
295
296 static std::shared_ptr<Grid> Shallow(std::shared_ptr<const Context> ctx, double Lx, double Ly,
297 double x0, double y0, unsigned int Mx, unsigned int My,
299
300 static std::shared_ptr<Grid> FromFile(std::shared_ptr<const Context> ctx,
301 const File &file,
302 const std::vector<std::string> &var_names,
304
305 static std::shared_ptr<Grid> FromOptions(std::shared_ptr<const Context> ctx);
306
307 std::shared_ptr<petsc::DM> get_dm(unsigned int dm_dof, unsigned int stencil_width) const;
308
309 std::shared_ptr<InputInterpolation> get_interpolation(const std::vector<double> &levels,
310 const File &input_file,
311 const std::string &variable_name,
312 InterpolationType type) const;
313
315
316 void report_parameters() const;
317
318 void compute_point_neighbors(double X, double Y,
319 int &i_left, int &i_right,
320 int &j_bottom, int &j_top) const;
321 std::vector<int> point_neighbors(double X, double Y) const;
322 std::vector<double> interpolation_weights(double x, double y) const;
323
324 unsigned int kBelowHeight(double height) const;
325
326 int max_patch_size() const;
327
328 std::shared_ptr<const Context> ctx() const;
329
330 int xs() const;
331 int xm() const;
332 int ys() const;
333 int ym() const;
334
335 const std::vector<double>& x() const;
336 double x(size_t i) const;
337
338 const std::vector<double>& y() const;
339 double y(size_t i) const;
340
341 const std::vector<double>& z() const;
342 double z(size_t i) const;
343
344 double dx() const;
345 double dy() const;
346 double cell_area() const;
347
348 unsigned int Mx() const;
349 unsigned int My() const;
350 unsigned int Mz() const;
351
352 double Lx() const;
353 double Ly() const;
354 double Lz() const;
355 double x0() const;
356 double y0() const;
357
358 const MappingInfo& get_mapping_info() const;
359 void set_mapping_info(const MappingInfo &info);
360
361 double dz_min() const;
362 double dz_max() const;
363
366
367 unsigned int size() const;
368 int rank() const;
369
370 const MPI_Comm com;
371
372 Vars& variables();
373 const Vars& variables() const;
374
375 PointsWithGhosts points(unsigned int stencil_width = 0) const {
376 return {*this, stencil_width};
377 }
378
379private:
380 struct Impl;
382
383 // Hide copy constructor / assignment operator.
384 Grid(const Grid &);
385 Grid & operator=(const Grid &);
386};
387
388namespace grid {
389
390std::vector<double> compute_vertical_levels(double new_Lz, size_t new_Mz,
391 grid::VerticalSpacing spacing, double Lambda = 0.0);
392
393//! Returns the distance from the point (i,j) to the origin.
394double radius(const Grid &grid, int i, int j);
395
396//! @brief Check if a point `(i,j)` is in the strip of `stripwidth` meters around the edge
397//! of the computational domain.
398inline bool in_null_strip(const Grid& grid, int i, int j, double strip_width) {
399 return (strip_width > 0.0 &&
400 (grid.x(i) <= grid.x(0) + strip_width ||
401 grid.x(i) >= grid.x(grid.Mx()-1) - strip_width ||
402 grid.y(j) <= grid.y(0) + strip_width ||
403 grid.y(j) >= grid.y(grid.My()-1) - strip_width));
404}
405
406/*!
407 * Return `true` if a point `(i, j)` is at an edge of the domain defined by the `grid`.
408 */
409inline bool domain_edge(const Grid &grid, int i, int j) {
410 return ((j == 0) or (j == (int)grid.My() - 1) or (i == 0) or (i == (int)grid.Mx() - 1));
411}
412
413std::array<unsigned, 2> nprocs(unsigned int size, unsigned int Mx,
414 unsigned int My);
415
416std::vector<unsigned int> ownership_ranges(unsigned int Mx, unsigned int Nx);
417
418} // namespace grid
419
420} // end of namespace pism
421
422#endif /* PISM_GRID_H */
A class for storing and accessing PISM configuration flags and parameters.
High-level PISM I/O class.
Definition File.hh:55
double y0() const
Y-coordinate of the center of the domain.
Definition Grid.cc:707
const std::vector< double > & x() const
X-coordinates.
Definition Grid.cc:623
const std::vector< double > & y() const
Y-coordinates.
Definition Grid.cc:633
unsigned int kBelowHeight(double height) const
Return the index k into zlevels[] so that zlevels[k] <= height < zlevels[k+1] and k < Mz.
Definition Grid.cc:275
unsigned int Mz() const
Number of vertical levels.
Definition Grid.cc:618
double cell_area() const
Definition Grid.cc:662
double Ly() const
Half-width of the computational domain.
Definition Grid.cc:692
int ys() const
Global starting index of this processor's subset.
Definition Grid.cc:593
grid::Periodicity periodicity() const
Return grid periodicity.
Definition Grid.cc:531
int max_patch_size() const
Definition Grid.cc:714
double Lx() const
Half-width of the computational domain.
Definition Grid.cc:687
PointsWithGhosts points(unsigned int stencil_width=0) const
Definition Grid.hh:375
double dx() const
Horizontal grid spacing.
Definition Grid.cc:653
const std::vector< double > & z() const
Z-coordinates within the ice.
Definition Grid.cc:643
int rank() const
MPI rank.
Definition Grid.cc:568
static std::shared_ptr< Grid > FromFile(std::shared_ptr< const Context > ctx, const File &file, const std::vector< std::string > &var_names, grid::Registration r)
Create a grid using one of variables in var_names in file.
Definition Grid.cc:250
static std::shared_ptr< Grid > Shallow(std::shared_ptr< const Context > ctx, double Lx, double Ly, double x0, double y0, unsigned int Mx, unsigned int My, grid::Registration r, grid::Periodicity p)
Initialize a uniform, shallow (3 z-levels) grid with half-widths (Lx,Ly) and Mx by My nodes.
Definition Grid.cc:134
Vars & variables()
Dictionary of variables (2D and 3D fields) associated with this grid.
Definition Grid.cc:578
std::shared_ptr< petsc::DM > get_dm(unsigned int dm_dof, unsigned int stencil_width) const
Get a PETSc DM ("distributed array manager") object for given dof (number of degrees of freedom per g...
Definition Grid.cc:510
unsigned int My() const
Total grid size in the Y direction.
Definition Grid.cc:613
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 Grid.cc:458
void report_parameters() const
Report grid parameters.
Definition Grid.cc:383
std::vector< double > interpolation_weights(double x, double y) const
Compute 4 interpolation weights necessary for linear interpolation from the current grid....
Definition Grid.cc:488
double x0() const
X-coordinate of the center of the domain.
Definition Grid.cc:702
grid::Registration registration() const
Definition Grid.cc:535
Impl * m_impl
Definition Grid.hh:381
std::shared_ptr< const Context > ctx() const
Return execution context this grid corresponds to.
Definition Grid.cc:540
unsigned int Mx() const
Total grid size in the X direction.
Definition Grid.cc:608
const MappingInfo & get_mapping_info() const
Definition Grid.cc:1577
unsigned int size() const
MPI communicator size.
Definition Grid.cc:573
std::shared_ptr< InputInterpolation > get_interpolation(const std::vector< double > &levels, const File &input_file, const std::string &variable_name, InterpolationType type) const
Definition Grid.cc:1586
std::vector< int > point_neighbors(double X, double Y) const
Definition Grid.cc:479
const MPI_Comm com
Definition Grid.hh:370
Grid & operator=(const Grid &)
double Lz() const
Height of the computational domain.
Definition Grid.cc:697
double dz_min() const
Minimum vertical spacing.
Definition Grid.cc:667
double dz_max() const
Maximum vertical spacing.
Definition Grid.cc:677
Grid(const Grid &)
double dy() const
Horizontal grid spacing.
Definition Grid.cc:658
int xs() const
Global starting index of this processor's subset.
Definition Grid.cc:588
void set_mapping_info(const MappingInfo &info)
Definition Grid.cc:1581
static std::shared_ptr< Grid > FromOptions(std::shared_ptr< const Context > ctx)
Create a grid using command-line options and (possibly) an input file.
Definition Grid.cc:1441
int xm() const
Width of this processor's sub-domain.
Definition Grid.cc:598
int ym() const
Width of this processor's sub-domain.
Definition Grid.cc:603
void forget_interpolations()
Definition Grid.cc:1605
Describes the PISM grid and the distribution of data across processors.
Definition Grid.hh:290
A basic logging class.
Definition Logger.hh:40
int j() const
Definition Grid.hh:193
int i() const
Definition Grid.hh:190
Points(const Grid &g)
Definition Grid.hh:227
A class for passing PISM variables from the core to other parts of the code (such as climate couplers...
Definition Vars.hh:42
std::vector< double > y
y coordinates
Definition Grid.hh:93
double Lx
domain half-width
Definition Grid.hh:82
double z_min
minimal value of the z dimension
Definition Grid.hh:86
unsigned int t_len
Definition Grid.hh:76
std::string filename
Definition Grid.hh:97
double Ly
domain half-height
Definition Grid.hh:84
double y0
y-coordinate of the domain center
Definition Grid.hh:80
double x0
x-coordinate of the domain center
Definition Grid.hh:78
std::vector< double > z
z coordinates
Definition Grid.hh:95
double z_max
maximal value of the z dimension
Definition Grid.hh:88
std::vector< double > x
x coordinates
Definition Grid.hh:91
std::map< std::string, AxisType > dimension_types
Definition Grid.hh:104
std::string variable_name
Definition Grid.hh:102
void report(const Logger &log, int threshold, std::shared_ptr< units::System > s) const
Definition Grid.cc:885
Contains parameters of an input file grid.
Definition Grid.hh:68
Periodicity periodicity
Grid periodicity.
Definition Grid.hh:163
std::vector< double > z
Vertical levels.
Definition Grid.hh:165
std::string variable_name
Name of the variable used to initialize the instance (empty if not used)
Definition Grid.hh:172
void vertical_grid_from_options(const Config &config)
Process -Mz and -Lz; set z;.
Definition Grid.cc:1326
double y0
Domain center in the Y direction.
Definition Grid.hh:155
double x0
Domain center in the X direction.
Definition Grid.hh:153
std::vector< unsigned int > procs_y
Processor ownership ranges in the Y direction.
Definition Grid.hh:169
void ownership_ranges_from_options(const Config &config, unsigned int size)
Re-compute ownership ranges. Uses current values of Mx and My.
Definition Grid.cc:1192
std::vector< unsigned int > procs_x
Processor ownership ranges in the X direction.
Definition Grid.hh:167
void horizontal_size_and_extent_from_options(const Config &config)
Process -Lx, -Ly, -x0, -y0; set Lx, Ly, x0, y0.
Definition Grid.cc:1298
Registration registration
Grid registration.
Definition Grid.hh:161
static Parameters FromGridDefinition(std::shared_ptr< units::System > unit_system, const File &file, const std::string &variable_name, Registration registration)
Definition Grid.cc:1050
void validate() const
Validate data members.
Definition Grid.cc:1335
unsigned int Mx
Number of grid points in the X direction.
Definition Grid.hh:157
double Ly
Domain half-width in the Y direction.
Definition Grid.hh:151
double Lx
Domain half-width in the X direction.
Definition Grid.hh:149
unsigned int My
Number of grid points in the Y direction.
Definition Grid.hh:159
Grid parameters; used to collect defaults before an Grid is allocated.
Definition Grid.hh:123
VerticalSpacing
Definition Grid.hh:53
@ EQUAL
Definition Grid.hh:53
@ UNKNOWN
Definition Grid.hh:53
@ QUADRATIC
Definition Grid.hh:53
std::string periodicity_to_string(Periodicity p)
Convert Periodicity to a STL string.
Definition Grid.cc:805
VerticalSpacing string_to_spacing(const std::string &keyword)
Convert an STL string to SpacingType.
Definition Grid.cc:820
std::vector< unsigned int > ownership_ranges(unsigned int Mx, unsigned int Nx)
Computes processor ownership ranges corresponding to equal area distribution among processors.
Definition Grid.cc:1426
Periodicity
Definition Grid.hh:54
@ XY_PERIODIC
Definition Grid.hh:54
@ Y_PERIODIC
Definition Grid.hh:54
@ X_PERIODIC
Definition Grid.hh:54
@ NOT_PERIODIC
Definition Grid.hh:54
Registration string_to_registration(const std::string &keyword)
Definition Grid.cc:844
Registration
Definition Grid.hh:56
@ CELL_CENTER
Definition Grid.hh:56
@ CELL_CORNER
Definition Grid.hh:56
Periodicity string_to_periodicity(const std::string &keyword)
Convert a string to Periodicity.
Definition Grid.cc:783
bool in_null_strip(const Grid &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 Grid.hh:398
bool domain_edge(const Grid &grid, int i, int j)
Definition Grid.hh:409
std::string spacing_to_string(VerticalSpacing s)
Convert SpacingType to an STL string.
Definition Grid.cc:834
std::vector< double > compute_vertical_levels(double new_Lz, size_t new_Mz, grid::VerticalSpacing spacing, double lambda)
Set the vertical levels in the ice according to values in Mz (number of levels), Lz (domain height),...
Definition Grid.cc:737
std::string registration_to_string(Registration registration)
Definition Grid.cc:857
std::array< unsigned, 2 > nprocs(unsigned int size, unsigned int Mx, unsigned int My)
Computes the number of processors in the X- and Y-directions.
Definition Grid.cc:1382
static const double g
Definition exactTestP.cc:36
std::shared_ptr< Grid > grid(std::shared_ptr< Context > ctx, char testname)
Definition pism.cc:176
Internal structures of Grid.
Definition Grid.cc:51
const double radius
Definition test_cube.c:18