tIceGrid.cc - 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.cc (45901B)
---
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 #include <cassert>
20
21 #include <map>
22 #include <numeric>
23 #include <petscsys.h>
24 #include <gsl/gsl_interp.h>
25
26 #include "IceGrid.hh"
27 #include "pism_utilities.hh"
28 #include "Time.hh"
29 #include "Time_Calendar.hh"
30 #include "ConfigInterface.hh"
31 #include "pism_options.hh"
32 #include "error_handling.hh"
33 #include "pism/util/io/File.hh"
34 #include "pism/util/Vars.hh"
35 #include "pism/util/Logger.hh"
36 #include "pism/util/projection.hh"
37 #include "pism/pism_config.hh"
38
39 #if (Pism_USE_PIO==1)
40 // Why do I need this???
41 #define _NETCDF
42 #include <pio.h>
43 #endif
44
45 namespace pism {
46
47 //! Internal structures of IceGrid.
48 struct IceGrid::Impl {
49 Impl(Context::ConstPtr ctx);
50
51 petsc::DM::Ptr create_dm(int da_dof, int stencil_width) const;
52 void set_ownership_ranges(const std::vector<unsigned int> &procs_x,
53 const std::vector<unsigned int> &procs_y);
54
55 void compute_horizontal_coordinates();
56
57 Context::ConstPtr ctx;
58
59 MappingInfo mapping_info;
60
61 // int to match types used by MPI
62 int rank;
63 int size;
64
65 //! @brief array containing lenghts (in the x-direction) of processor sub-domains
66 std::vector<PetscInt> procs_x;
67 //! @brief array containing lenghts (in the y-direction) of processor sub-domains
68 std::vector<PetscInt> procs_y;
69
70 Periodicity periodicity;
71
72 GridRegistration registration;
73
74 //! x-coordinates of grid points
75 std::vector<double> x;
76 //! y-coordinates of grid points
77 std::vector<double> y;
78 //! vertical grid levels in the ice; correspond to the storage grid
79 std::vector<double> z;
80
81 int xs, xm, ys, ym;
82 //! horizontal grid spacing
83 double dx;
84 //! horizontal grid spacing
85 double dy;
86 //! cell area (meters^2)
87 double cell_area;
88 //! number of grid points in the x-direction
89 unsigned int Mx;
90 //! number of grid points in the y-direction
91 unsigned int My;
92
93 //! x-coordinate of the grid center
94 double x0;
95 //! y-coordinate of the grid center
96 double y0;
97
98 //! half width of the ice model grid in x-direction (m)
99 double Lx;
100 //! half width of the ice model grid in y-direction (m)
101 double Ly;
102
103 std::map<int,petsc::DM::WeakPtr> dms;
104
105 // This DM is used for I/O operations and is not owned by any
106 // IceModelVec (so far, anyway). We keep a pointer to it here to
107 // avoid re-allocating it many times.
108 petsc::DM::Ptr dm_scalar_global;
109
110 //! @brief A dictionary with pointers to IceModelVecs, for passing
111 //! them from the one component to another (e.g. from IceModel to
112 //! surface and ocean models).
113 Vars variables;
114
115 //! GSL binary search accelerator used to speed up kBelowHeight().
116 gsl_interp_accel *bsearch_accel;
117
118 //! ParallelIO I/O decompositions.
119 std::map<int, int> io_decompositions;
120 };
121
122 IceGrid::Impl::Impl(Context::ConstPtr context)
123 : ctx(context), mapping_info("mapping", ctx->unit_system()) {
124 // empty
125 }
126
127 //! Convert a string to Periodicity.
128 Periodicity string_to_periodicity(const std::string &keyword) {
129 if (keyword == "none") {
130 return NOT_PERIODIC;
131 } else if (keyword == "x") {
132 return X_PERIODIC;
133 } else if (keyword == "y") {
134 return Y_PERIODIC;
135 } else if (keyword == "xy") {
136 return XY_PERIODIC;
137 } else {
138 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "grid periodicity type '%s' is invalid.",
139 keyword.c_str());
140 }
141 }
142
143 //! Convert Periodicity to a STL string.
144 std::string periodicity_to_string(Periodicity p) {
145 switch (p) {
146 case NOT_PERIODIC:
147 return "none";
148 case X_PERIODIC:
149 return "x";
150 case Y_PERIODIC:
151 return "y";
152 default:
153 case XY_PERIODIC:
154 return "xy";
155 }
156 }
157
158 //! Convert an STL string to SpacingType.
159 SpacingType string_to_spacing(const std::string &keyword) {
160 if (keyword == "quadratic") {
161 return QUADRATIC;
162 } else if (keyword == "equal") {
163 return EQUAL;
164 } else {
165 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "ice vertical spacing type '%s' is invalid.",
166 keyword.c_str());
167 }
168 }
169
170 //! Convert SpacingType to an STL string.
171 std::string spacing_to_string(SpacingType s) {
172 switch (s) {
173 case EQUAL:
174 return "equal";
175 default:
176 case QUADRATIC:
177 return "quadratic";
178 }
179 }
180
181 GridRegistration string_to_registration(const std::string &keyword) {
182 if (keyword == "center") {
183 return CELL_CENTER;
184 } else if (keyword == "corner") {
185 return CELL_CORNER;
186 } else {
187 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid grid registration: %s",
188 keyword.c_str());
189 }
190 }
191
192 std::string registration_to_string(GridRegistration registration) {
193 switch (registration) {
194 case CELL_CORNER:
195 return "corner";
196 default:
197 case CELL_CENTER:
198 return "center";
199 }
200 }
201
202 /*! @brief Initialize a uniform, shallow (3 z-levels) grid with half-widths (Lx,Ly) and Mx by My
203 * nodes.
204 */
205 IceGrid::Ptr IceGrid::Shallow(Context::ConstPtr ctx,
206 double Lx, double Ly,
207 double x0, double y0,
208 unsigned int Mx, unsigned int My,
209 GridRegistration registration,
210 Periodicity periodicity) {
211 try {
212 GridParameters p(ctx->config());
213 p.Lx = Lx;
214 p.Ly = Ly;
215 p.x0 = x0;
216 p.y0 = y0;
217 p.Mx = Mx;
218 p.My = My;
219 p.registration = registration;
220 p.periodicity = periodicity;
221
222 double Lz = ctx->config()->get_number("grid.Lz");
223 p.z.resize(3);
224 p.z[0] = 0.0;
225 p.z[1] = 0.5 * Lz;
226 p.z[2] = 1.0 * Lz;
227
228 p.ownership_ranges_from_options(ctx->size());
229
230 return IceGrid::Ptr(new IceGrid(ctx, p));
231 } catch (RuntimeError &e) {
232 e.add_context("initializing a shallow grid");
233 throw;
234 }
235 }
236
237 //! @brief Create a PISM distributed computational grid.
238 IceGrid::IceGrid(Context::ConstPtr context, const GridParameters &p)
239 : com(context->com()), m_impl(new Impl(context)) {
240
241 try {
242 m_impl->bsearch_accel = gsl_interp_accel_alloc();
243 if (m_impl->bsearch_accel == NULL) {
244 throw RuntimeError(PISM_ERROR_LOCATION, "Failed to allocate a GSL interpolation accelerator");
245 }
246
247 MPI_Comm_rank(com, &m_impl->rank);
248 MPI_Comm_size(com, &m_impl->size);
249
250 p.validate();
251
252 m_impl->Mx = p.Mx;
253 m_impl->My = p.My;
254 m_impl->x0 = p.x0;
255 m_impl->y0 = p.y0;
256 m_impl->Lx = p.Lx;
257 m_impl->Ly = p.Ly;
258 m_impl->registration = p.registration;
259 m_impl->periodicity = p.periodicity;
260 m_impl->z = p.z;
261 m_impl->set_ownership_ranges(p.procs_x, p.procs_y);
262
263 m_impl->compute_horizontal_coordinates();
264
265 {
266 unsigned int stencil_width = (unsigned int)context->config()->get_number("grid.max_stencil_width");
267
268 try {
269 petsc::DM::Ptr tmp = this->get_dm(1, stencil_width);
270 } catch (RuntimeError &e) {
271 e.add_context("distributing a %d x %d grid across %d processors.",
272 Mx(), My(), size());
273 throw;
274 }
275
276 // hold on to a DM corresponding to dof=1, stencil_width=0 (it will
277 // be needed for I/O operations)
278 m_impl->dm_scalar_global = this->get_dm(1, 0);
279
280 DMDALocalInfo info;
281 PetscErrorCode ierr = DMDAGetLocalInfo(*m_impl->dm_scalar_global, &info);
282 PISM_CHK(ierr, "DMDAGetLocalInfo");
283
284 m_impl->xs = info.xs;
285 m_impl->xm = info.xm;
286 m_impl->ys = info.ys;
287 m_impl->ym = info.ym;
288
289 }
290 } catch (RuntimeError &e) {
291 e.add_context("allocating IceGrid");
292 throw;
293 }
294 }
295
296 //! Create a grid using one of variables in `var_names` in `file`.
297 IceGrid::Ptr IceGrid::FromFile(Context::ConstPtr ctx,
298 const std::string &filename,
299 const std::vector<std::string> &var_names,
300 GridRegistration r) {
301
302 File file(ctx->com(), filename, PISM_NETCDF3, PISM_READONLY);
303
304 for (auto name : var_names) {
305 if (file.find_variable(name)) {
306 return FromFile(ctx, file, name, r);
307 }
308 }
309
310 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "file %s does not have any of %s."
311 " Cannot initialize the grid.",
312 filename.c_str(),
313 join(var_names, ",").c_str());
314 }
315
316 //! Create a grid from a file, get information from variable `var_name`.
317 IceGrid::Ptr IceGrid::FromFile(Context::ConstPtr ctx,
318 const File &file,
319 const std::string &var_name,
320 GridRegistration r) {
321 try {
322 const Logger &log = *ctx->log();
323
324 // The following call may fail because var_name does not exist. (And this is fatal!)
325 // Note that this sets defaults using configuration parameters, too.
326 GridParameters p(ctx, file, var_name, r);
327
328 // if we have no vertical grid information, create a fake 2-level vertical grid.
329 if (p.z.size() < 2) {
330 double Lz = ctx->config()->get_number("grid.Lz");
331 log.message(3,
332 "WARNING: Can't determine vertical grid information using '%s' in %s'\n"
333 " Using 2 levels and Lz of %3.3fm\n",
334 var_name.c_str(), file.filename().c_str(), Lz);
335
336 p.z = {0.0, Lz};
337 }
338
339
340 p.ownership_ranges_from_options(ctx->size());
341
342 return IceGrid::Ptr(new IceGrid(ctx, p));
343 } catch (RuntimeError &e) {
344 e.add_context("initializing computational grid from variable \"%s\" in \"%s\"",
345 var_name.c_str(), file.filename().c_str());
346 throw;
347 }
348 }
349
350 IceGrid::~IceGrid() {
351 gsl_interp_accel_free(m_impl->bsearch_accel);
352
353 #if (Pism_USE_PIO==1)
354 for (auto p : m_impl->io_decompositions) {
355 int ierr = PIOc_freedecomp(m_impl->ctx->pio_iosys_id(), p.second);
356 if (ierr != PIO_NOERR) {
357 m_impl->ctx->log()->message(1, "Failed to de-allocate a ParallelIO decomposition");
358 }
359 }
360 #endif
361
362 delete m_impl;
363 }
364
365 //! \brief Set the vertical levels in the ice according to values in `Mz` (number of levels), `Lz`
366 //! (domain height), `spacing` (quadratic or equal) and `lambda` (quadratic spacing parameter).
367 /*!
368 - When `vertical_spacing == EQUAL`, the vertical grid in the ice is equally spaced:
369 `zlevels[k] = k dz` where `dz = Lz / (Mz - 1)`.
370 - When `vertical_spacing == QUADRATIC`, the spacing is a quadratic function. The intent
371 is that the spacing is smaller near the base than near the top.
372
373 In particular, if
374 \f$\zeta_k = k / (\mathtt{Mz} - 1)\f$ then `zlevels[k] = Lz *`
375 ((\f$\zeta_k\f$ / \f$\lambda\f$) * (1.0 + (\f$\lambda\f$ - 1.0)
376 * \f$\zeta_k\f$)) where \f$\lambda\f$ = 4. The value \f$\lambda\f$
377 indicates the slope of the quadratic function as it leaves the base.
378 Thus a value of \f$\lambda\f$ = 4 makes the spacing about four times finer
379 at the base than equal spacing would be.
380 */
381 std::vector<double> IceGrid::compute_vertical_levels(double new_Lz, unsigned int new_Mz,
382 SpacingType spacing, double lambda) {
383
384 if (new_Mz < 2) {
385 throw RuntimeError(PISM_ERROR_LOCATION, "Mz must be at least 2");
386 }
387
388 if (new_Lz <= 0) {
389 throw RuntimeError(PISM_ERROR_LOCATION, "Lz must be positive");
390 }
391
392 if (spacing == QUADRATIC and lambda <= 0) {
393 throw RuntimeError(PISM_ERROR_LOCATION, "lambda must be positive");
394 }
395
396 std::vector<double> result(new_Mz);
397
398 // Fill the levels in the ice:
399 switch (spacing) {
400 case EQUAL: {
401 double dz = new_Lz / ((double) new_Mz - 1);
402
403 // Equal spacing
404 for (unsigned int k=0; k < new_Mz - 1; k++) {
405 result[k] = dz * ((double) k);
406 }
407 result[new_Mz - 1] = new_Lz; // make sure it is exactly equal
408 break;
409 }
410 case QUADRATIC: {
411 // this quadratic scheme is an attempt to be less extreme in the fineness near the base.
412 for (unsigned int k=0; k < new_Mz - 1; k++) {
413 const double zeta = ((double) k) / ((double) new_Mz - 1);
414 result[k] = new_Lz * ((zeta / lambda) * (1.0 + (lambda - 1.0) * zeta));
415 }
416 result[new_Mz - 1] = new_Lz; // make sure it is exactly equal
417 break;
418 }
419 default:
420 throw RuntimeError(PISM_ERROR_LOCATION, "spacing can not be UNKNOWN");
421 }
422
423 return result;
424 }
425
426 //! Return the index `k` into `zlevels[]` so that `zlevels[k] <= height < zlevels[k+1]` and `k < Mz`.
427 unsigned int IceGrid::kBelowHeight(double height) const {
428
429 if (height < 0.0 - 1.0e-6) {
430 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "height = %5.4f is below base of ice"
431 " (height must be non-negative)\n", height);
432 }
433
434 if (height > Lz() + 1.0e-6) {
435 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "height = %5.4f is above top of computational"
436 " grid Lz = %5.4f\n", height, Lz());
437 }
438
439 return gsl_interp_accel_find(m_impl->bsearch_accel, &m_impl->z[0], m_impl->z.size(), height);
440 }
441
442 //! \brief Computes the number of processors in the X- and Y-directions.
443 static void compute_nprocs(unsigned int Mx, unsigned int My, unsigned int size,
444 unsigned int &Nx, unsigned int &Ny) {
445
446 if (My <= 0) {
447 throw RuntimeError(PISM_ERROR_LOCATION, "'My' is invalid.");
448 }
449
450 Nx = (unsigned int)(0.5 + sqrt(((double)Mx)*((double)size)/((double)My)));
451 Ny = 0;
452
453 if (Nx == 0) {
454 Nx = 1;
455 }
456
457 while (Nx > 0) {
458 Ny = size / Nx;
459 if (Nx*Ny == (unsigned int)size) {
460 break;
461 }
462 Nx--;
463 }
464
465 if (Mx > My and Nx < Ny) {
466 // Swap Nx and Ny
467 int tmp = Nx;
468 Nx = Ny;
469 Ny = tmp;
470 }
471
472 if ((Mx / Nx) < 2) { // note: integer division
473 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't split %d grid points into %d parts (X-direction).",
474 Mx, (int)Nx);
475 }
476
477 if ((My / Ny) < 2) { // note: integer division
478 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't split %d grid points into %d parts (Y-direction).",
479 My, (int)Ny);
480 }
481 }
482
483
484 //! \brief Computes processor ownership ranges corresponding to equal area
485 //! distribution among processors.
486 static std::vector<unsigned int> ownership_ranges(unsigned int Mx,
487 unsigned int Nx) {
488
489 std::vector<unsigned int> result(Nx);
490
491 for (unsigned int i=0; i < Nx; i++) {
492 result[i] = Mx / Nx + ((Mx % Nx) > i);
493 }
494 return result;
495 }
496
497 //! Set processor ownership ranges. Takes care of type conversion (`unsigned int` -> `PetscInt`).
498 void IceGrid::Impl::set_ownership_ranges(const std::vector<unsigned int> &input_procs_x,
499 const std::vector<unsigned int> &input_procs_y) {
500 if (input_procs_x.size() * input_procs_y.size() != (size_t)size) {
501 throw RuntimeError(PISM_ERROR_LOCATION, "length(procs_x) * length(procs_y) != MPI size");
502 }
503
504 procs_x.resize(input_procs_x.size());
505 for (unsigned int k = 0; k < input_procs_x.size(); ++k) {
506 procs_x[k] = input_procs_x[k];
507 }
508
509 procs_y.resize(input_procs_y.size());
510 for (unsigned int k = 0; k < input_procs_y.size(); ++k) {
511 procs_y[k] = input_procs_y[k];
512 }
513 }
514
515 struct OwnershipRanges {
516 std::vector<unsigned int> x, y;
517 };
518
519 //! Compute processor ownership ranges using the grid size, MPI communicator size, and command-line
520 //! options `-Nx`, `-Ny`, `-procs_x`, `-procs_y`.
521 static OwnershipRanges compute_ownership_ranges(unsigned int Mx,
522 unsigned int My,
523 unsigned int size) {
524 OwnershipRanges result;
525
526 unsigned int Nx_default, Ny_default;
527 compute_nprocs(Mx, My, size, Nx_default, Ny_default);
528
529 // check -Nx and -Ny
530 options::Integer Nx("-Nx", "Number of processors in the x direction", Nx_default);
531 options::Integer Ny("-Ny", "Number of processors in the y direction", Ny_default);
532
533 // validate results (compute_nprocs checks its results, but we also need to validate command-line
534 // options)
535 if ((Mx / Nx) < 2) {
536 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't split %d grid points into %d parts (X-direction).",
537 Mx, (int)Nx);
538 }
539
540 if ((My / Ny) < 2) {
541 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't split %d grid points into %d parts (Y-direction).",
542 My, (int)Ny);
543 }
544
545 if (Nx * Ny != (int)size) {
546 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Nx * Ny has to be equal to %d.", size);
547 }
548
549 // check -procs_x and -procs_y
550 options::IntegerList procs_x("-procs_x", "Processor ownership ranges (x direction)", {});
551 options::IntegerList procs_y("-procs_y", "Processor ownership ranges (y direction)", {});
552
553 if (procs_x.is_set()) {
554 if (procs_x->size() != (unsigned int)Nx) {
555 throw RuntimeError(PISM_ERROR_LOCATION, "-Nx has to be equal to the -procs_x size.");
556 }
557
558 result.x.resize(procs_x->size());
559 for (unsigned int k = 0; k < procs_x->size(); ++k) {
560 result.x[k] = procs_x[k];
561 }
562
563 } else {
564 result.x = ownership_ranges(Mx, Nx);
565 }
566
567 if (procs_y.is_set()) {
568 if (procs_y->size() != (unsigned int)Ny) {
569 throw RuntimeError(PISM_ERROR_LOCATION, "-Ny has to be equal to the -procs_y size.");
570 }
571
572 result.y.resize(procs_y->size());
573 for (unsigned int k = 0; k < procs_y->size(); ++k) {
574 result.y[k] = procs_y[k];
575 }
576 } else {
577 result.y = ownership_ranges(My, Ny);
578 }
579
580 if (result.x.size() * result.y.size() != size) {
581 throw RuntimeError(PISM_ERROR_LOCATION, "length(procs_x) * length(procs_y) != MPI size");
582 }
583
584 return result;
585 }
586
587 //! Compute horizontal grid spacing. See compute_horizontal_coordinates() for more.
588 static double compute_horizontal_spacing(double half_width, unsigned int M, bool cell_centered) {
589 if (cell_centered) {
590 return 2.0 * half_width / M;
591 } else {
592 return 2.0 * half_width / (M - 1);
593 }
594 }
595
596 //! Compute grid coordinates for one direction (X or Y).
597 static std::vector<double> compute_coordinates(unsigned int M, double delta,
598 double v_min, double v_max,
599 bool cell_centered) {
600 std::vector<double> result(M);
601
602 // Here v_min, v_max define the extent of the computational domain,
603 // which is not necessarily the same thing as the smallest and
604 // largest values of grid coordinates.
605 if (cell_centered) {
606 for (unsigned int i = 0; i < M; ++i) {
607 result[i] = v_min + (i + 0.5) * delta;
608 }
609 result[M - 1] = v_max - 0.5*delta;
610 } else {
611 for (unsigned int i = 0; i < M; ++i) {
612 result[i] = v_min + i * delta;
613 }
614 result[M - 1] = v_max;
615 }
616 return result;
617 }
618
619 //! Compute horizontal spacing parameters `dx` and `dy` and grid coordinates using `Mx`, `My`, `Lx`, `Ly` and periodicity.
620 /*!
621 The grid used in PISM, in particular the PETSc DAs used here, are periodic in x and y.
622 This means that the ghosted values ` foo[i+1][j], foo[i-1][j], foo[i][j+1], foo[i][j-1]`
623 for all 2D Vecs, and similarly in the x and y directions for 3D Vecs, are always available.
624 That is, they are available even if i,j is a point at the edge of the grid. On the other
625 hand, by default, `dx` is the full width `2 * Lx` divided by `Mx - 1`.
626 This means that we conceive of the computational domain as starting at the `i = 0`
627 grid location and ending at the `i = Mx - 1` grid location, in particular.
628 This idea is not quite compatible with the periodic nature of the grid.
629
630 The upshot is that if one computes in a truly periodic way then the gap between the
631 `i = 0` and `i = Mx - 1` grid points should \em also have width `dx`.
632 Thus we compute `dx = 2 * Lx / Mx`.
633 */
634 void IceGrid::Impl::compute_horizontal_coordinates() {
635
636 bool cell_centered = registration == CELL_CENTER;
637
638 dx = compute_horizontal_spacing(Lx, Mx, cell_centered);
639
640 dy = compute_horizontal_spacing(Ly, My, cell_centered);
641
642 cell_area = dx * dy;
643
644 double
645 x_min = x0 - Lx,
646 x_max = x0 + Lx;
647
648 x = compute_coordinates(Mx, dx,
649 x_min, x_max,
650 cell_centered);
651
652 double
653 y_min = y0 - Ly,
654 y_max = y0 + Ly;
655
656 y = compute_coordinates(My, dy,
657 y_min, y_max,
658 cell_centered);
659 }
660
661 //! \brief Report grid parameters.
662 void IceGrid::report_parameters() const {
663
664 const Logger &log = *this->ctx()->log();
665 units::System::Ptr sys = this->ctx()->unit_system();
666
667 log.message(2, "computational domain and grid:\n");
668
669 units::Converter km(sys, "m", "km");
670
671 // report on grid
672 log.message(2,
673 " grid size %d x %d x %d\n",
674 Mx(), My(), Mz());
675
676 // report on computational box
677 log.message(2,
678 " spatial domain %.2f km x %.2f km x %.2f m\n",
679 km(2*Lx()), km(2*Ly()), Lz());
680
681 // report on grid cell dims
682 if ((dx() && dy()) > 1000.) {
683 log.message(2,
684 " horizontal grid cell %.2f km x %.2f km\n",
685 km(dx()), km(dy()));
686 } else {
687 log.message(2,
688 " horizontal grid cell %.0f m x %.0f m\n",
689 dx(), dy());
690 }
691 if (fabs(dz_max() - dz_min()) <= 1.0e-8) {
692 log.message(2,
693 " vertical spacing in ice dz = %.3f m (equal spacing)\n",
694 dz_min());
695 } else {
696 log.message(2,
697 " vertical spacing in ice uneven, %d levels, %.3f m < dz < %.3f m\n",
698 Mz(), dz_min(), dz_max());
699 }
700
701 // if -verbose (=-verbose 3) then (somewhat redundantly) list parameters of grid
702 {
703 log.message(3,
704 " IceGrid parameters:\n");
705 log.message(3,
706 " Lx = %6.2f km, Ly = %6.2f km, Lz = %6.2f m, \n",
707 km(Lx()), km(Ly()), Lz());
708 log.message(3,
709 " x0 = %6.2f km, y0 = %6.2f km, (coordinates of center)\n",
710 km(x0()), km(y0()));
711 log.message(3,
712 " Mx = %d, My = %d, Mz = %d, \n",
713 Mx(), My(), Mz());
714 log.message(3,
715 " dx = %6.3f km, dy = %6.3f km, \n",
716 km(dx()), km(dy()));
717 log.message(3,
718 " Nx = %d, Ny = %d\n",
719 (int)m_impl->procs_x.size(), (int)m_impl->procs_y.size());
720
721 log.message(3,
722 " Registration: %s\n",
723 registration_to_string(m_impl->registration).c_str());
724 log.message(3,
725 " Periodicity: %s\n",
726 periodicity_to_string(m_impl->periodicity).c_str());
727 }
728
729 {
730 log.message(5,
731 " REALLY verbose output on IceGrid:\n");
732 log.message(5,
733 " vertical levels in ice (Mz=%d, Lz=%5.4f): ", Mz(), Lz());
734 for (unsigned int k=0; k < Mz(); k++) {
735 log.message(5, " %5.4f, ", z(k));
736 }
737 log.message(5, "\n");
738 }
739
740 }
741
742
743 //! \brief Computes indices of grid points to the lower left and upper right from (X,Y).
744 /*!
745 * \code
746 * 3 2
747 * o-------o
748 * | |
749 * | + |
750 * o-------o
751 * 0 1
752 * \endcode
753 *
754 * If "+" is the point (X,Y), then (i_left, j_bottom) corresponds to
755 * point "0" and (i_right, j_top) corresponds to point "2".
756 *
757 * Does not check if the resulting indexes are in the current
758 * processor's domain. Ensures that computed indexes are within the
759 * grid.
760 */
761 void IceGrid::compute_point_neighbors(double X, double Y,
762 int &i_left, int &i_right,
763 int &j_bottom, int &j_top) const {
764 i_left = (int)floor((X - m_impl->x[0])/m_impl->dx);
765 j_bottom = (int)floor((Y - m_impl->y[0])/m_impl->dy);
766
767 i_right = i_left + 1;
768 j_top = j_bottom + 1;
769
770 i_left = std::max(i_left, 0);
771 i_right = std::max(i_right, 0);
772
773 i_left = std::min(i_left, (int)m_impl->Mx - 1);
774 i_right = std::min(i_right, (int)m_impl->Mx - 1);
775
776 j_bottom = std::max(j_bottom, 0);
777 j_top = std::max(j_top, 0);
778
779 j_bottom = std::min(j_bottom, (int)m_impl->My - 1);
780 j_top = std::min(j_top, (int)m_impl->My - 1);
781 }
782
783 std::vector<int> IceGrid::compute_point_neighbors(double X, double Y) const {
784 int i_left, i_right, j_bottom, j_top;
785 this->compute_point_neighbors(X, Y, i_left, i_right, j_bottom, j_top);
786 return {i_left, i_right, j_bottom, j_top};
787 }
788
789 //! \brief Compute 4 interpolation weights necessary for linear interpolation
790 //! from the current grid. See compute_point_neighbors for the ordering of
791 //! neighbors.
792 std::vector<double> IceGrid::compute_interp_weights(double X, double Y) const{
793 int i_left = 0, i_right = 0, j_bottom = 0, j_top = 0;
794 // these values (zeros) are used when interpolation is impossible
795 double alpha = 0.0, beta = 0.0;
796
797 compute_point_neighbors(X, Y, i_left, i_right, j_bottom, j_top);
798
799 if (i_left != i_right) {
800 assert(m_impl->x[i_right] - m_impl->x[i_left] != 0.0);
801 alpha = (X - m_impl->x[i_left]) / (m_impl->x[i_right] - m_impl->x[i_left]);
802 }
803
804 if (j_bottom != j_top) {
805 assert(m_impl->y[j_top] - m_impl->y[j_bottom] != 0.0);
806 beta = (Y - m_impl->y[j_bottom]) / (m_impl->y[j_top] - m_impl->y[j_bottom]);
807 }
808
809 return {(1.0 - alpha) * (1.0 - beta),
810 alpha * (1.0 - beta),
811 alpha * beta,
812 (1.0 - alpha) * beta};
813 }
814
815 // Computes the hash corresponding to the DM with given dof and stencil_width.
816 static int dm_hash(int dm_dof, int stencil_width) {
817 if (dm_dof < 0 or dm_dof > IceGrid::max_dm_dof) {
818 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
819 "Invalid dm_dof argument: %d", dm_dof);
820 }
821
822 if (stencil_width < 0 or stencil_width > IceGrid::max_stencil_width) {
823 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
824 "Invalid stencil_width argument: %d", stencil_width);
825 }
826
827 return IceGrid::max_stencil_width * dm_dof + stencil_width;
828 }
829
830 //! @brief Get a PETSc DM ("distributed array manager") object for given `dof` (number of degrees of
831 //! freedom per grid point) and stencil width.
832 petsc::DM::Ptr IceGrid::get_dm(int da_dof, int stencil_width) const {
833 petsc::DM::Ptr result;
834
835 int j = dm_hash(da_dof, stencil_width);
836
837 if (m_impl->dms[j].expired()) {
838 result = m_impl->create_dm(da_dof, stencil_width);
839 m_impl->dms[j] = result;
840 } else {
841 result = m_impl->dms[j].lock();
842 }
843
844 return result;
845 }
846
847 //! Return grid periodicity.
848 Periodicity IceGrid::periodicity() const {
849 return m_impl->periodicity;
850 }
851
852 GridRegistration IceGrid::registration() const {
853 return m_impl->registration;
854 }
855
856 //! Return execution context this grid corresponds to.
857 Context::ConstPtr IceGrid::ctx() const {
858 return m_impl->ctx;
859 }
860
861 //! @brief Create a DM with the given number of `dof` (degrees of freedom per grid point) and
862 //! stencil width.
863 petsc::DM::Ptr IceGrid::Impl::create_dm(int da_dof, int stencil_width) const {
864
865 ctx->log()->message(3,
866 "* Creating a DM with dof=%d and stencil_width=%d...\n",
867 da_dof, stencil_width);
868
869 DM result;
870 PetscErrorCode ierr = DMDACreate2d(ctx->com(),
871 DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,
872 DMDA_STENCIL_BOX,
873 Mx, My,
874 (PetscInt)procs_x.size(),
875 (PetscInt)procs_y.size(),
876 da_dof, stencil_width,
877 &procs_x[0], &procs_y[0], // lx, ly
878 &result);
879 PISM_CHK(ierr,"DMDACreate2d");
880
881 #if PETSC_VERSION_GE(3,8,0)
882 ierr = DMSetUp(result); PISM_CHK(ierr,"DMSetUp");
883 #endif
884
885 return petsc::DM::Ptr(new petsc::DM(result));
886 }
887
888 //! MPI rank.
889 int IceGrid::rank() const {
890 return m_impl->rank;
891 }
892
893 //! MPI communicator size.
894 unsigned int IceGrid::size() const {
895 return m_impl->size;
896 }
897
898 //! Dictionary of variables (2D and 3D fields) associated with this grid.
899 Vars& IceGrid::variables() {
900 return m_impl->variables;
901 }
902
903 //! Dictionary of variables (2D and 3D fields) associated with this grid.
904 const Vars& IceGrid::variables() const {
905 return m_impl->variables;
906 }
907
908 //! Global starting index of this processor's subset.
909 int IceGrid::xs() const {
910 return m_impl->xs;
911 }
912
913 //! Global starting index of this processor's subset.
914 int IceGrid::ys() const {
915 return m_impl->ys;
916 }
917
918 //! Width of this processor's sub-domain.
919 int IceGrid::xm() const {
920 return m_impl->xm;
921 }
922
923 //! Width of this processor's sub-domain.
924 int IceGrid::ym() const {
925 return m_impl->ym;
926 }
927
928 //! Total grid size in the X direction.
929 unsigned int IceGrid::Mx() const {
930 return m_impl->Mx;
931 }
932
933 //! Total grid size in the Y direction.
934 unsigned int IceGrid::My() const {
935 return m_impl->My;
936 }
937
938 //! Number of vertical levels.
939 unsigned int IceGrid::Mz() const {
940 return m_impl->z.size();
941 }
942
943 //! X-coordinates.
944 const std::vector<double>& IceGrid::x() const {
945 return m_impl->x;
946 }
947
948 //! Get a particular x coordinate.
949 double IceGrid::x(size_t i) const {
950 return m_impl->x[i];
951 }
952
953 //! Y-coordinates.
954 const std::vector<double>& IceGrid::y() const {
955 return m_impl->y;
956 }
957
958 //! Get a particular y coordinate.
959 double IceGrid::y(size_t i) const {
960 return m_impl->y[i];
961 }
962
963 //! Z-coordinates within the ice.
964 const std::vector<double>& IceGrid::z() const {
965 return m_impl->z;
966 }
967
968 //! Get a particular z coordinate.
969 double IceGrid::z(size_t i) const {
970 return m_impl->z[i];
971 }
972
973 //! Horizontal grid spacing.
974 double IceGrid::dx() const {
975 return m_impl->dx;
976 }
977
978 //! Horizontal grid spacing.
979 double IceGrid::dy() const {
980 return m_impl->dy;
981 }
982
983 double IceGrid::cell_area() const {
984 return m_impl->cell_area;
985 }
986
987 //! Minimum vertical spacing.
988 double IceGrid::dz_min() const {
989 double result = m_impl->z.back();
990 for (unsigned int k = 0; k < m_impl->z.size() - 1; ++k) {
991 const double dz = m_impl->z[k + 1] - m_impl->z[k];
992 result = std::min(dz, result);
993 }
994 return result;
995 }
996
997 //! Maximum vertical spacing.
998 double IceGrid::dz_max() const {
999 double result = 0.0;
1000 for (unsigned int k = 0; k < m_impl->z.size() - 1; ++k) {
1001 const double dz = m_impl->z[k + 1] - m_impl->z[k];
1002 result = std::max(dz, result);
1003 }
1004 return result;
1005 }
1006
1007 //! Half-width of the computational domain.
1008 double IceGrid::Lx() const {
1009 return m_impl->Lx;
1010 }
1011
1012 //! Half-width of the computational domain.
1013 double IceGrid::Ly() const {
1014 return m_impl->Ly;
1015 }
1016
1017 //! Height of the computational domain.
1018 double IceGrid::Lz() const {
1019 return m_impl->z.back();
1020 }
1021
1022 //! X-coordinate of the center of the domain.
1023 double IceGrid::x0() const {
1024 return m_impl->x0;
1025 }
1026
1027 //! Y-coordinate of the center of the domain.
1028 double IceGrid::y0() const {
1029 return m_impl->y0;
1030 }
1031
1032 //! @brief Returns the distance from the point (i,j) to the origin.
1033 double radius(const IceGrid &grid, int i, int j) {
1034 return sqrt(grid.x(i) * grid.x(i) + grid.y(j) * grid.y(j));
1035 }
1036
1037 // grid_info
1038
1039 void grid_info::reset() {
1040
1041 filename = "";
1042
1043 t_len = 0;
1044 time = 0;
1045
1046 x_len = 0;
1047 x0 = 0;
1048 Lx = 0;
1049
1050 y_len = 0;
1051 y0 = 0;
1052 Ly = 0;
1053
1054 z_len = 0;
1055 z_min = 0;
1056 z_max = 0;
1057 }
1058
1059 grid_info::grid_info() {
1060 reset();
1061 }
1062
1063 void grid_info::report(const Logger &log, int threshold, units::System::Ptr s) const {
1064 units::Converter km(s, "m", "km");
1065
1066 log.message(threshold,
1067 " x: %5d points, [%10.3f, %10.3f] km, x0 = %10.3f km, Lx = %10.3f km\n",
1068 this->x_len,
1069 km(this->x0 - this->Lx),
1070 km(this->x0 + this->Lx),
1071 km(this->x0),
1072 km(this->Lx));
1073
1074 log.message(threshold,
1075 " y: %5d points, [%10.3f, %10.3f] km, y0 = %10.3f km, Ly = %10.3f km\n",
1076 this->y_len,
1077 km(this->y0 - this->Ly),
1078 km(this->y0 + this->Ly),
1079 km(this->y0),
1080 km(this->Ly));
1081
1082 log.message(threshold,
1083 " z: %5d points, [%10.3f, %10.3f] m\n",
1084 this->z_len, this->z_min, this->z_max);
1085
1086 log.message(threshold,
1087 " t: %5d points, last time = %.3f years\n\n",
1088 this->t_len, units::convert(s, this->time, "seconds", "years"));
1089 }
1090
1091 grid_info::grid_info(const File &file, const std::string &variable,
1092 units::System::Ptr unit_system,
1093 GridRegistration r) {
1094 try {
1095 reset();
1096
1097 filename = file.filename();
1098
1099 // try "variable" as the standard_name first, then as the short name:
1100 auto var = file.find_variable(variable, variable);
1101
1102 if (not var.exists) {
1103 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "variable \"%s\" is missing", variable.c_str());
1104 }
1105
1106 auto dimensions = file.dimensions(var.name);
1107
1108 for (auto dimension_name : dimensions) {
1109
1110 AxisType dimtype = file.dimension_type(dimension_name, unit_system);
1111
1112 switch (dimtype) {
1113 case X_AXIS:
1114 {
1115 this->x = file.read_dimension(dimension_name);
1116 this->x_len = this->x.size();
1117 double
1118 x_min = vector_min(this->x),
1119 x_max = vector_max(this->x);
1120 this->x0 = 0.5 * (x_min + x_max);
1121 this->Lx = 0.5 * (x_max - x_min);
1122 if (r == CELL_CENTER) {
1123 const double dx = this->x[1] - this->x[0];
1124 this->Lx += 0.5 * dx;
1125 }
1126 break;
1127 }
1128 case Y_AXIS:
1129 {
1130 this->y = file.read_dimension(dimension_name);
1131 this->y_len = this->y.size();
1132 double
1133 y_min = vector_min(this->y),
1134 y_max = vector_max(this->y);
1135 this->y0 = 0.5 * (y_min + y_max);
1136 this->Ly = 0.5 * (y_max - y_min);
1137 if (r == CELL_CENTER) {
1138 const double dy = this->y[1] - this->y[0];
1139 this->Ly += 0.5 * dy;
1140 }
1141 break;
1142 }
1143 case Z_AXIS:
1144 {
1145 this->z = file.read_dimension(dimension_name);
1146 this->z_len = this->z.size();
1147 this->z_min = vector_min(this->z);
1148 this->z_max = vector_max(this->z);
1149 break;
1150 }
1151 case T_AXIS:
1152 {
1153 this->t_len = file.dimension_length(dimension_name);
1154 this->time = vector_max(file.read_dimension(dimension_name));
1155 break;
1156 }
1157 default:
1158 {
1159 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
1160 "can't figure out which direction dimension '%s' corresponds to.",
1161 dimension_name.c_str());
1162 }
1163 } // switch
1164 } // for loop
1165 } catch (RuntimeError &e) {
1166 e.add_context("getting grid information using variable '%s' in '%s'", variable.c_str(),
1167 file.filename().c_str());
1168 throw;
1169 }
1170 }
1171
1172 GridParameters::GridParameters() {
1173
1174 // set to something invalid
1175 Lx = -1.0;
1176 Ly = -1.0;
1177
1178 x0 = 0.0;
1179 y0 = 0.0;
1180
1181 Mx = 0;
1182 My = 0;
1183
1184 registration = CELL_CENTER;
1185 periodicity = NOT_PERIODIC;
1186 }
1187
1188 GridParameters::GridParameters(Config::ConstPtr config) {
1189 init_from_config(config);
1190 }
1191
1192 void GridParameters::ownership_ranges_from_options(unsigned int size) {
1193 OwnershipRanges procs = compute_ownership_ranges(Mx, My, size);
1194 procs_x = procs.x;
1195 procs_y = procs.y;
1196 }
1197
1198 //! Initialize from a configuration database. Does not try to compute ownership ranges.
1199 void GridParameters::init_from_config(Config::ConstPtr config) {
1200 Lx = config->get_number("grid.Lx");
1201 Ly = config->get_number("grid.Ly");
1202
1203 x0 = 0.0;
1204 y0 = 0.0;
1205
1206 Mx = config->get_number("grid.Mx");
1207 My = config->get_number("grid.My");
1208
1209 periodicity = string_to_periodicity(config->get_string("grid.periodicity"));
1210 registration = string_to_registration(config->get_string("grid.registration"));
1211
1212 double Lz = config->get_number("grid.Lz");
1213 unsigned int Mz = config->get_number("grid.Mz");
1214 double lambda = config->get_number("grid.lambda");
1215 SpacingType s = string_to_spacing(config->get_string("grid.ice_vertical_spacing"));
1216 z = IceGrid::compute_vertical_levels(Lz, Mz, s, lambda);
1217 // does not set ownership ranges because we don't know if these settings are final
1218 }
1219
1220 void GridParameters::init_from_file(Context::ConstPtr ctx,
1221 const File &file,
1222 const std::string &variable_name,
1223 GridRegistration r) {
1224 int size = 0;
1225 MPI_Comm_size(ctx->com(), &size);
1226
1227 // set defaults (except for ownership ranges) from configuration parameters
1228 init_from_config(ctx->config());
1229
1230 grid_info input_grid(file, variable_name, ctx->unit_system(), r);
1231
1232 Lx = input_grid.Lx;
1233 Ly = input_grid.Ly;
1234 x0 = input_grid.x0;
1235 y0 = input_grid.y0;
1236 Mx = input_grid.x_len;
1237 My = input_grid.y_len;
1238 registration = r;
1239 z = input_grid.z;
1240 }
1241
1242 GridParameters::GridParameters(Context::ConstPtr ctx,
1243 const File &file,
1244 const std::string &variable_name,
1245 GridRegistration r) {
1246 init_from_file(ctx, file, variable_name, r);
1247 }
1248
1249 GridParameters::GridParameters(Context::ConstPtr ctx,
1250 const std::string &filename,
1251 const std::string &variable_name,
1252 GridRegistration r) {
1253 File file(ctx->com(), filename, PISM_NETCDF3, PISM_READONLY);
1254 init_from_file(ctx, file, variable_name, r);
1255 }
1256
1257
1258 void GridParameters::horizontal_size_from_options() {
1259 Mx = options::Integer("-Mx", "grid size in X direction", Mx);
1260 My = options::Integer("-My", "grid size in Y direction", My);
1261 }
1262
1263 void GridParameters::horizontal_extent_from_options() {
1264 // Domain size
1265 {
1266 Lx = 1000.0 * options::Real("-Lx", "Half of the grid extent in the Y direction, in km",
1267 Lx / 1000.0);
1268 Ly = 1000.0 * options::Real("-Ly", "Half of the grid extent in the X direction, in km",
1269 Ly / 1000.0);
1270 }
1271
1272 // Alternatively: domain size and extent
1273 {
1274 options::RealList x_range("-x_range", "min,max x coordinate values", {});
1275 options::RealList y_range("-y_range", "min,max y coordinate values", {});
1276
1277 if (x_range.is_set() and y_range.is_set()) {
1278 if (x_range->size() != 2 or y_range->size() != 2) {
1279 throw RuntimeError(PISM_ERROR_LOCATION, "-x_range and/or -y_range argument is invalid.");
1280 }
1281 x0 = (x_range[0] + x_range[1]) / 2.0;
1282 y0 = (y_range[0] + y_range[1]) / 2.0;
1283 Lx = (x_range[1] - x_range[0]) / 2.0;
1284 Ly = (y_range[1] - y_range[0]) / 2.0;
1285 }
1286 }
1287 }
1288
1289 void GridParameters::vertical_grid_from_options(Config::ConstPtr config) {
1290 double Lz = z.size() > 0 ? z.back() : config->get_number("grid.Lz");
1291 int Mz = z.size() > 0 ? z.size() : config->get_number("grid.Mz");
1292 double lambda = config->get_number("grid.lambda");
1293 SpacingType s = string_to_spacing(config->get_string("grid.ice_vertical_spacing"));
1294
1295 z = IceGrid::compute_vertical_levels(Lz, Mz, s, lambda);
1296 }
1297
1298 void GridParameters::validate() const {
1299 if (Mx < 3) {
1300 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Mx = %d is invalid (has to be 3 or greater)", Mx);
1301 }
1302
1303 if (My < 3) {
1304 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "My = %d is invalid (has to be 3 or greater)", My);
1305 }
1306
1307 if (Lx <= 0.0) {
1308 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Lx = %f is invalid (negative)", Lx);
1309 }
1310
1311 if (Ly <= 0.0) {
1312 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Ly = %f is invalid (negative)", Ly);
1313 }
1314
1315 if (not is_increasing(z)) {
1316 throw RuntimeError(PISM_ERROR_LOCATION, "z levels are not increasing");
1317 }
1318
1319 if (z[0] > 1e-6) {
1320 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "first z level is not zero: %f", z[0]);
1321 }
1322
1323 if (z.back() < 0.0) {
1324 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "last z level is negative: %f", z.back());
1325 }
1326
1327 if (std::accumulate(procs_x.begin(), procs_x.end(), 0.0) != Mx) {
1328 throw RuntimeError(PISM_ERROR_LOCATION, "procs_x don't sum up to Mx");
1329 }
1330
1331 if (std::accumulate(procs_y.begin(), procs_y.end(), 0.0) != My) {
1332 throw RuntimeError(PISM_ERROR_LOCATION, "procs_y don't sum up to My");
1333 }
1334 }
1335
1336 //! Create a grid using command-line options and (possibly) an input file.
1337 /** Processes options -i, -bootstrap, -Mx, -My, -Mz, -Lx, -Ly, -Lz, -x_range, -y_range.
1338 */
1339 IceGrid::Ptr IceGrid::FromOptions(Context::ConstPtr ctx) {
1340 auto config = ctx->config();
1341
1342 auto input_file = config->get_string("input.file");
1343 bool bootstrap = config->get_flag("input.bootstrap");
1344
1345 GridRegistration r = string_to_registration(config->get_string("grid.registration"));
1346
1347 Logger::ConstPtr log = ctx->log();
1348
1349 if (not input_file.empty() and (not bootstrap)) {
1350 // These options are ignored because we're getting *all* the grid
1351 // parameters from a file.
1352 options::ignored(*log, "-Mx");
1353 options::ignored(*log, "-My");
1354 options::ignored(*log, "-Mz");
1355 options::ignored(*log, "-Mbz");
1356 options::ignored(*log, "-Lx");
1357 options::ignored(*log, "-Ly");
1358 options::ignored(*log, "-Lz");
1359 options::ignored(*log, "-z_spacing");
1360
1361 // get grid from a PISM input file
1362 return IceGrid::FromFile(ctx, input_file, {"enthalpy", "temp"}, r);
1363 } else if (not input_file.empty() and bootstrap) {
1364 // bootstrapping; get domain size defaults from an input file, allow overriding all grid
1365 // parameters using command-line options
1366
1367 GridParameters input_grid(config);
1368
1369 std::vector<std::string> names = {"land_ice_thickness", "bedrock_altitude",
1370 "thk", "topg"};
1371 bool grid_info_found = false;
1372
1373 File file(ctx->com(), input_file, PISM_NETCDF3, PISM_READONLY);
1374
1375 for (auto name : names) {
1376
1377 grid_info_found = file.find_variable(name);
1378 if (not grid_info_found) {
1379 // Failed to find using a short name. Try using name as a
1380 // standard name...
1381 grid_info_found = file.find_variable("unlikely_name", name).exists;
1382 }
1383
1384 if (grid_info_found) {
1385 input_grid = GridParameters(ctx, file, name, r);
1386 break;
1387 }
1388 }
1389
1390 if (not grid_info_found) {
1391 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no geometry information found in '%s'",
1392 input_file.c_str());
1393 }
1394
1395 // process all possible options controlling grid parameters, overriding values read from a file
1396 input_grid.horizontal_size_from_options();
1397 input_grid.horizontal_extent_from_options();
1398 input_grid.vertical_grid_from_options(config);
1399 input_grid.ownership_ranges_from_options(ctx->size());
1400
1401 IceGrid::Ptr result(new IceGrid(ctx, input_grid));
1402
1403 units::System::Ptr sys = ctx->unit_system();
1404 units::Converter km(sys, "m", "km");
1405
1406 // report on resulting computational box
1407 log->message(2,
1408 " setting computational box for ice from '%s' and\n"
1409 " user options: [%6.2f km, %6.2f km] x [%6.2f km, %6.2f km] x [0 m, %6.2f m]\n",
1410 input_file.c_str(),
1411 km(result->x0() - result->Lx()),
1412 km(result->x0() + result->Lx()),
1413 km(result->y0() - result->Ly()),
1414 km(result->y0() + result->Ly()),
1415 result->Lz());
1416
1417 return result;
1418 } else {
1419 // This covers the two remaining cases "-i is not set, -bootstrap is set" and "-i is not set,
1420 // -bootstrap is not set either".
1421 throw RuntimeError(PISM_ERROR_LOCATION,
1422 "Please set the input file using the \"-i\" command-line option.");
1423 }
1424 }
1425
1426 const MappingInfo& IceGrid::get_mapping_info() const {
1427 return m_impl->mapping_info;
1428 }
1429
1430 void IceGrid::set_mapping_info(const MappingInfo &info) {
1431 m_impl->mapping_info = info;
1432 // FIXME: re-compute lat/lon coordinates
1433 }
1434
1435 #if (Pism_USE_PIO==1)
1436 static int pio_decomp_hash(int dof, int output_datatype) {
1437 if (dof < 0 or dof > IceGrid::max_dm_dof) {
1438 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
1439 "Invalid dof argument: %d", dof);
1440 }
1441
1442 // pio.h includes netcdf.h and netcdf.h defines NC_FIRSTUSERTYPEID which exceeds
1443 // all constants corresponding to built-in types
1444 return NC_FIRSTUSERTYPEID * dof + output_datatype;
1445 }
1446 #endif
1447
1448 /*!
1449 * initialize an I/O decomposition
1450 *
1451 * @param[in] dof size of the last dimension (usually z)
1452 * @param[in] output_datatype an integer specifying a data type (`PIO_DOUBLE`, etc)
1453 */
1454 int IceGrid::pio_io_decomposition(int dof, int output_datatype) const {
1455 int result = 0;
1456 #if (Pism_USE_PIO==1)
1457 {
1458 int hash = pio_decomp_hash(dof, output_datatype);
1459 result = m_impl->io_decompositions[hash];
1460
1461 if (result == 0) {
1462
1463 int ndims = dof < 2 ? 2 : 3;
1464
1465 // the last element is not used if ndims == 2
1466 std::vector<int> gdimlen{(int)My(), (int)Mx(), dof};
1467 std::vector<long int> start{ys(), xs(), 0}, count{ym(), xm(), dof};
1468
1469 int stat = PIOc_InitDecomp_bc(m_impl->ctx->pio_iosys_id(),
1470 output_datatype, ndims, gdimlen.data(),
1471 start.data(), count.data(), &result);
1472 m_impl->io_decompositions[hash] = result;
1473 if (stat != PIO_NOERR) {
1474 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
1475 "Failed to create a ParallelIO I/O decomposition");
1476 }
1477 }
1478 }
1479 #else
1480 (void) dof;
1481 (void) output_datatype;
1482 #endif
1483 return result;
1484 }
1485
1486 } // end of namespace pism