tIceGrid_Regional.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_Regional.cc (5911B)
---
1 /* Copyright (C) 2015, 2016, 2017, 2018, 2019 PISM Authors
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
20 #include <vector>
21 #include <cassert>
22
23 #include <gsl/gsl_interp.h>
24
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/IceGrid.hh"
28 #include "pism/util/io/File.hh"
29 #include "pism/util/Component.hh" // process_input_options
30
31 namespace pism {
32
33 static void validate_range(const std::string &axis,
34 const std::vector<double> &x,
35 double x_min, double x_max) {
36 if (x_min >= x_max) {
37 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid -%s_range: %s_min >= %s_max (%f >= %f).",
38 axis.c_str(), axis.c_str(), axis.c_str(),
39 x_min, x_max);
40 }
41
42 if (x_min >= x.back()) {
43 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid -%s_range: %s_min >= max(%s) (%f >= %f).",
44 axis.c_str(), axis.c_str(), axis.c_str(),
45 x_min, x.back());
46 }
47
48 if (x_max <= x.front()) {
49 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid -%s-range: %s_max <= min(%s) (%f <= %f).",
50 axis.c_str(), axis.c_str(), axis.c_str(),
51 x_max, x.front());
52 }
53 }
54
55 static void subset_extent(const std::string& axis,
56 const std::vector<double> &x,
57 double x_min, double x_max,
58 double &x0, double &Lx, unsigned int &Mx) {
59
60 validate_range(axis, x, x_min, x_max);
61
62 size_t x_start = gsl_interp_bsearch(&x[0], x_min, 0, x.size() - 1);
63 // include one more point if we can
64 if (x_start > 0) {
65 x_start -= 1;
66 }
67
68 size_t x_end = gsl_interp_bsearch(&x[0], x_max, 0, x.size() - 1);
69 // include one more point if we can
70 x_end = std::min(x.size() - 1, x_end + 1);
71
72 // NOTE: this assumes the CELL_CORNER grid registration
73 Lx = (x[x_end] - x[x_start]) / 2.0;
74
75 x0 = (x[x_start] + x[x_end]) / 2.0;
76
77 assert(x_end + 1 >= x_start);
78 Mx = x_end + 1 - x_start;
79 }
80
81 //! Create a grid using command-line options and (possibly) an input file.
82 /** Processes options -i, -bootstrap, -Mx, -My, -Mz, -Lx, -Ly, -Lz, -x_range, -y_range.
83 */
84 IceGrid::Ptr regional_grid_from_options(Context::Ptr ctx) {
85
86 auto options = process_input_options(ctx->com(), ctx->config());
87
88 const options::RealList x_range("-x_range",
89 "range of X coordinates in the selected subset", {});
90 const options::RealList y_range("-y_range",
91 "range of Y coordinates in the selected subset", {});
92
93 if (options.type == INIT_BOOTSTRAP and x_range.is_set() and y_range.is_set()) {
94 // bootstrapping; get domain size defaults from an input file, allow overriding all grid
95 // parameters using command-line options
96
97 if (x_range->size() != 2) {
98 throw RuntimeError(PISM_ERROR_LOCATION, "invalid -x_range argument: need 2 numbers.");
99 }
100
101 if (y_range->size() != 2) {
102 throw RuntimeError(PISM_ERROR_LOCATION, "invalid -y_range argument: need 2 numbers.");
103 }
104
105 GridParameters input_grid(ctx->config());
106
107 std::vector<std::string> names = {"enthalpy", "temp", "land_ice_thickness",
108 "bedrock_altitude", "thk", "topg"};
109 bool grid_info_found = false;
110
111 File file(ctx->com(), options.filename, PISM_NETCDF3, PISM_READONLY);
112 for (auto name : names) {
113
114 grid_info_found = file.find_variable(name);
115 if (not grid_info_found) {
116 // Failed to find using a short name. Try using name as a
117 // standard name...
118 grid_info_found = file.find_variable("unlikely_name", name).exists;
119 }
120
121 if (grid_info_found) {
122 input_grid = GridParameters(ctx, file, name, CELL_CORNER);
123
124 grid_info full = grid_info(file, name, ctx->unit_system(), CELL_CORNER);
125
126 // x direction
127 subset_extent("x", full.x, x_range[0], x_range[1],
128 input_grid.x0, input_grid.Lx, input_grid.Mx);
129 // y direction
130 subset_extent("y", full.y, y_range[0], y_range[1],
131 input_grid.y0, input_grid.Ly, input_grid.My);
132
133 // Set registration to "CELL_CORNER" so that IceGrid computes
134 // coordinates correctly.
135 input_grid.registration = CELL_CORNER;
136
137 break;
138 }
139 }
140
141 if (not grid_info_found) {
142 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
143 "no geometry information found in '%s'",
144 options.filename.c_str());
145 }
146
147 // process options controlling vertical grid parameters, overriding values read from a file
148 input_grid.vertical_grid_from_options(ctx->config());
149
150 // process options controlling ownership ranges
151 input_grid.ownership_ranges_from_options(ctx->size());
152
153 return IceGrid::Ptr(new IceGrid(ctx, input_grid));
154 } else if (x_range.is_set() ^ y_range.is_set()) {
155 throw RuntimeError(PISM_ERROR_LOCATION, "Please set both -x_range and -y_range.");
156 } else {
157 return IceGrid::FromOptions(ctx);
158 }
159 }
160
161
162 } // end of namespace pism