Loading [MathJax]/jax/output/HTML-CSS/config.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
Array.hh
Go to the documentation of this file.
1// Copyright (C) 2008--2023 Ed Bueler, Constantine Khroulev, and 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#ifndef PISM_ARRAY_H
20#define PISM_ARRAY_H
21
22#include <initializer_list>
23#include <memory> // shared_ptr, dynamic_pointer_cast
24#include <cstdint> // uint64_t
25
26#include "pism/util/error_handling.hh" // RuntimeError
27
28namespace pism {
29
31
32namespace io {
33enum Type : int;
34class Default;
35}
36
37class Grid;
38class File;
40
41namespace petsc {
42class DM;
43class Vec;
44class Viewer;
45} // end of namespace petsc
46
47namespace units {
48class System;
49}
50
52public:
53 virtual ~PetscAccessible() = default;
54 virtual void begin_access() const = 0;
55 virtual void end_access() const = 0;
56};
57
58namespace array {
59
60//! What "kind" of a vector to create: with or without ghosts.
62
63//! Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
65public:
67 AccessScope(std::initializer_list<const PetscAccessible *> vecs);
70 void add(const PetscAccessible &v);
71 void add(const std::vector<const PetscAccessible*> &vecs);
72private:
73 std::vector<const PetscAccessible*> m_vecs;
74};
75
76/*!
77 * Interpolation helper. Does not check if points needed for interpolation are within the current
78 * processor's sub-domain.
79 */
80template<class F, typename T>
81T interpolate(const F &field, double x, double y) {
82 auto grid = field.grid();
83
84 int i_left = 0, i_right = 0, j_bottom = 0, j_top = 0;
85 grid->compute_point_neighbors(x, y, i_left, i_right, j_bottom, j_top);
86
87 auto w = grid->compute_interp_weights(x, y);
88
89 return (w[0] * field(i_left, j_bottom) +
90 w[1] * field(i_right, j_bottom) +
91 w[2] * field(i_right, j_top) +
92 w[3] * field(i_left, j_top));
93}
94
95//! \brief Abstract class for reading, writing, allocating, and accessing a
96//! DA-based PETSc Vec (2D and 3D fields) from within IceModel.
97/*!
98 @anchor icemodelvec_use
99
100 This class represents 2D and 3D fields in PISM. Its methods common to all
101 the derived classes can be split (roughly) into six kinds:
102
103 - memory allocation (create)
104 - point-wise access (begin_access(), end_access())
105 - arithmetic (range(), norm(), add(), shift(), scale(), set(), ...)
106 - setting or reading metadata (set_attrs(), metadata())
107 - file input/output (read, write, regrid)
108 - tracking whether a field was updated (get_state_counter(), inc_state_counter())
109
110 ## Memory allocation
111
112 \code
113 array::Scalar var(grid, "var_name", WITH_GHOSTS);
114 // var is ready to use
115 \endcode
116
117 ("WITH_GHOSTS" means "can be used in computations using map-plane neighbors
118 of grid points.)
119
120 It is usually a good idea to set variable metadata right after creating it.
121 The method set_attrs() is used throughout PISM to set commonly used
122 attributes.
123
124 ## Point-wise access
125
126 PETSc performs some pointer arithmetic magic to allow convenient indexing of
127 grid point values. Because of this one needs to surround the code using row,
128 column or level indexes with begin_access() and end_access() calls:
129
130 \code
131 double foo;
132 int i = 0, j = 0;
133 array::Scalar var;
134 // assume that var was allocated
135 ierr = var.begin_access(); CHKERRQ(ierr);
136 foo = var(i,j) * 2;
137 ierr = var.end_access(); CHKERRQ(ierr);
138 \endcode
139
140 Please see [this page](@ref computational_grid) for a discussion of the
141 organization of PISM's computational grid and examples of for-loops you will
142 probably put between begin_access() and end_access().
143
144 To ensure that ghost values are up to date add the following call
145 before the code using ghosts:
146
147 \code
148 ierr = var.update_ghosts(); CHKERRQ(ierr);
149 \endcode
150
151 ## Reading and writing variables
152
153 PISM can read variables either from files with data on a grid matching the
154 current grid (read()) or, using bilinear interpolation, from files
155 containing data on a different (but compatible) grid (regrid()).
156
157 To write a field to a "prepared" NetCDF file, use write(). (A file is prepared
158 if it contains all the necessary dimensions, coordinate variables and global
159 metadata.)
160
161 If you need to "prepare" a file, do:
162 \code
163 File file(grid.com, PISM_NETCDF3);
164 io::prepare_for_output(file, *grid.ctx());
165 \endcode
166
167 A note about NetCDF write performance: due to limitations of the NetCDF
168 (classic, version 3) format, it is significantly faster to
169 \code
170 for (all variables)
171 var.define(...);
172
173 for (all variables)
174 var.write(...);
175 \endcode
176
177 as opposed to
178
179 \code
180 for (all variables) {
181 var.define(...);
182 var.write(...);
183 }
184 \endcode
185
186 array::Array::define() is here so that we can use the first approach.
187
188 ## Tracking if a field changed
189
190 It is possible to track if a certain field changed with the help of
191 get_state_counter() and inc_state_counter() methods.
192
193 For example, PISM's SIA code re-computes the smoothed bed only if the bed
194 deformation code updated it:
195
196 \code
197 if (bed->get_state_counter() > bed_state_counter) {
198 ierr = bed_smoother->preprocess_bed(...); CHKERRQ(ierr);
199 bed_state_counter = bed->get_state_counter();
200 }
201 \endcode
202
203 The state counter is **not** updated automatically. For the code snippet above
204 to work, a bed deformation model has to call inc_state_counter() after an
205 update.
206*/
207class Array : public PetscAccessible {
208public:
209 virtual ~Array();
210
211 std::shared_ptr<const Grid> grid() const;
212 unsigned int ndims() const;
213 std::vector<int> shape() const;
214 //! @brief Returns the number of degrees of freedom per grid point.
215 unsigned int ndof() const;
216 unsigned int stencil_width() const;
217
218 const std::vector<double>& levels() const;
219
220 std::vector<double> norm(int n) const;
221
222 void add(double alpha, const Array &x);
223 void shift(double alpha);
224 void scale(double alpha);
225
226 petsc::Vec& vec() const;
227 std::shared_ptr<petsc::DM> dm() const;
228
229 void set_name(const std::string &name);
230 const std::string& get_name() const;
231
232 void define(const File &file, io::Type default_type) const;
233
234 void read(const std::string &filename, unsigned int time);
235 void read(const File &file, unsigned int time);
236
237 void write(const std::string &filename) const;
238 void write(const File &file) const;
239
240 void regrid(const std::string &filename, io::Default default_value);
241 void regrid(const File &file, io::Default default_value);
242
243 virtual void begin_access() const;
244 virtual void end_access() const;
245 void update_ghosts();
246
247 std::shared_ptr<petsc::Vec> allocate_proc0_copy() const;
248 void put_on_proc0(petsc::Vec &onp0) const;
249 void get_from_proc0(petsc::Vec &onp0);
250
251 void set(double c);
252
253 SpatialVariableMetadata& metadata(unsigned int N = 0);
254
255 const SpatialVariableMetadata& metadata(unsigned int N = 0) const;
256
257 int state_counter() const;
258 void inc_state_counter();
259
261
262 void view(std::vector<std::shared_ptr<petsc::Viewer> > viewers) const;
263
264protected:
265 Array(std::shared_ptr<const Grid> grid,
266 const std::string &name,
267 Kind ghostedp,
268 size_t dof,
269 size_t stencil_width,
270 const std::vector<double> &zlevels);
271 struct Impl;
273
274 // will be cast to double**, double***, or Vector2** in derived classes
275 // This is not hidden in m_impl to make it possible to inline operator()
276 mutable void *m_array;
277
278 void set_begin_access_use_dof(bool flag);
279
280 void read_impl(const File &file, unsigned int time);
281 virtual void regrid_impl(const File &file, io::Default default_value);
282 void write_impl(const File &file) const;
283
284 void checkCompatibility(const char *function, const Array &other) const;
285
286 //! @brief Check array indices and warn if they are out of range.
287 void check_array_indices(int i, int j, unsigned int k) const;
288
289 void copy_to_vec(std::shared_ptr<petsc::DM> destination_da, petsc::Vec &destination) const;
290 void get_dof(std::shared_ptr<petsc::DM> da_result, petsc::Vec &result, unsigned int start,
291 unsigned int count=1) const;
292 void set_dof(std::shared_ptr<petsc::DM> da_source, petsc::Vec &source, unsigned int start,
293 unsigned int count=1);
294private:
295 size_t size() const;
296 // disable copy constructor and the assignment operator:
297 Array(const Array &other);
299public:
300 //! Dump an Array to a file. *This is for debugging only.*
301 //! Uses const char[] to make it easier to call it from gdb.
302 void dump(const char filename[]) const;
303
304 uint64_t fletcher64_serial() const;
305 uint64_t fletcher64() const;
306 std::string checksum(bool serial) const;
307 void print_checksum(const char *prefix = "", bool serial=false) const;
308
309protected:
310 void put_on_proc0(petsc::Vec &parallel, petsc::Vec &onp0) const;
311 void get_from_proc0(petsc::Vec &onp0, petsc::Vec &parallel) const;
312};
313
314//! `std::dynamic_pointer_cast` wrapper that checks if the cast succeeded.
315template <class T>
316static typename std::shared_ptr<T> cast(std::shared_ptr<Array> input) {
317 auto result = std::dynamic_pointer_cast<T, Array>(input);
318
319 if (not result) {
320 throw RuntimeError(PISM_ERROR_LOCATION, "dynamic cast failure");
321 }
322
323 return result;
324}
325
326} // end of namespace array
327
328/**
329 * Convert a PETSc Vec from the units in `from` into units in `to` (in place).
330 *
331 * @param v data to convert
332 * @param system unit system
333 * @param spec1 source unit specification string
334 * @param spec2 destination unit specification string
335 */
336void convert_vec(petsc::Vec &v, std::shared_ptr<units::System> system,
337 const std::string &spec1, const std::string &spec2);
338
339} // end of namespace pism
340
341#endif /* PISM_ARRAY_H */
High-level PISM I/O class.
Definition File.hh:55
Describes the PISM grid and the distribution of data across processors.
Definition Grid.hh:290
virtual ~PetscAccessible()=default
virtual void begin_access() const =0
virtual void end_access() const =0
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
void add(const PetscAccessible &v)
Definition Array.cc:896
std::vector< const PetscAccessible * > m_vecs
Definition Array.hh:73
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void read(const std::string &filename, unsigned int time)
Definition Array.cc:731
void write_impl(const File &file) const
Writes an Array to a NetCDF file.
Definition Array.cc:487
virtual void end_access() const
Checks if an Array is allocated and calls DAVecRestoreArray.
Definition Array.cc:589
virtual void regrid_impl(const File &file, io::Default default_value)
Gets an Array from a file file, interpolating onto the current grid.
Definition Array.cc:378
unsigned int ndof() const
Returns the number of degrees of freedom per grid point.
Definition Array.cc:135
void set_interpolation_type(InterpolationType type)
Definition Array.cc:178
void set_name(const std::string &name)
Sets the variable name to name.
Definition Array.cc:342
void dump(const char filename[]) const
Dumps a variable to a file, overwriting this file's contents (for debugging).
Definition Array.cc:530
Array(const Array &other)
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition Array.cc:463
virtual void begin_access() const
Checks if an Array is allocated and calls DAVecGetArray.
Definition Array.cc:568
void shift(double alpha)
Result: v[j] <- v[j] + alpha for all j. Calls VecShift.
Definition Array.cc:216
const std::vector< double > & levels() const
Definition Array.cc:139
size_t size() const
Return the total number of elements in the owned part of an array.
Definition Array.cc:909
void get_from_proc0(petsc::Vec &onp0)
Gets a local Array2 from processor 0.
Definition Array.cc:1058
void write(const std::string &filename) const
Definition Array.cc:722
void copy_to_vec(std::shared_ptr< petsc::DM > destination_da, petsc::Vec &destination) const
Copies v to a global vector 'destination'. Ghost points are discarded.
Definition Array.cc:238
const std::string & get_name() const
Get the name of an Array object.
Definition Array.cc:354
virtual ~Array()
Definition Array.cc:104
petsc::Vec & vec() const
Definition Array.cc:310
uint64_t fletcher64_serial() const
Definition Array.cc:1075
void get_dof(std::shared_ptr< petsc::DM > da_result, petsc::Vec &result, unsigned int start, unsigned int count=1) const
Definition Array.cc:248
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition Array.cc:224
Array & operator=(const Array &)
void print_checksum(const char *prefix="", bool serial=false) const
Definition Array.cc:1161
int state_counter() const
Get the object state counter.
Definition Array.cc:127
std::vector< int > shape() const
Definition Array.cc:159
std::string checksum(bool serial) const
Definition Array.cc:1152
void view(std::vector< std::shared_ptr< petsc::Viewer > > viewers) const
View a 2D vector field using existing PETSc viewers.
Definition Array.cc:1168
std::shared_ptr< const Grid > grid() const
Definition Array.cc:131
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void set_dof(std::shared_ptr< petsc::DM > da_source, petsc::Vec &source, unsigned int start, unsigned int count=1)
Definition Array.cc:273
void inc_state_counter()
Increment the object state counter.
Definition Array.cc:150
std::shared_ptr< petsc::DM > dm() const
Definition Array.cc:324
unsigned int ndims() const
Returns the number of spatial dimensions.
Definition Array.cc:155
void regrid(const std::string &filename, io::Default default_value)
Definition Array.cc:736
std::shared_ptr< petsc::Vec > allocate_proc0_copy() const
Definition Array.cc:926
uint64_t fletcher64() const
Definition Array.cc:1107
void put_on_proc0(petsc::Vec &onp0) const
Puts a local array::Scalar on processor 0.
Definition Array.cc:1015
void set_begin_access_use_dof(bool flag)
Definition Array.cc:174
void add(double alpha, const Array &x)
Result: v <- v + alpha * x. Calls VecAXPY.
Definition Array.cc:206
std::vector< double > norm(int n) const
Computes the norm of all the components of an Array.
Definition Array.cc:668
void update_ghosts()
Updates ghost points.
Definition Array.cc:615
unsigned int stencil_width() const
Get the stencil width of the current Array. Returns 0 if ghosts are not available.
Definition Array.cc:302
void checkCompatibility(const char *function, const Array &other) const
Checks if two Arrays have compatible sizes, dimensions and numbers of degrees of freedom.
Definition Array.cc:545
void read_impl(const File &file, unsigned int time)
Reads appropriate NetCDF variable(s) into an Array.
Definition Array.cc:416
void check_array_indices(int i, int j, unsigned int k) const
Check array indices and warn if they are out of range.
Definition Array.cc:636
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Definition Array.hh:207
#define PISM_ERROR_LOCATION
#define n
Definition exactTestM.c:37
static std::shared_ptr< T > cast(std::shared_ptr< Array > input)
std::dynamic_pointer_cast wrapper that checks if the cast succeeded.
Definition Array.hh:316
Kind
What "kind" of a vector to create: with or without ghosts.
Definition Array.hh:61
@ WITH_GHOSTS
Definition Array.hh:61
@ WITHOUT_GHOSTS
Definition Array.hh:61
T interpolate(const F &field, double x, double y)
Definition Array.hh:81
static double F(double SL, double B, double H, double alpha)
static const double k
Definition exactTestP.cc:42
void convert_vec(petsc::Vec &v, units::System::Ptr system, const std::string &spec1, const std::string &spec2)
Definition Array.cc:1225
int count
Definition test_cube.c:16