tgrid.rst - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tgrid.rst (10636B)
---
1 .. include:: ../../../global.txt
2
3 .. contents::
4
5 .. _sec-grid:
6
7 Spatial grid
8 ------------
9
10 The PISM grid covering the computational box is equally spaced in horizontal (`x`
11 and `y`) directions. Vertical spacing in the ice is quadratic by default but
12 optionally equal spacing can be chosen; choose with options :opt:`-z_spacing`
13 [``quadratic``, ``equal``\] at bootstrapping. The grid read from a "``-i``" input file is
14 used as is. The bedrock thermal layer model always uses equal vertical spacing.
15
16 The grid is described by four numbers, namely the number of grid points ``Mx`` in the
17 `x` direction, the number ``My`` in the `y` direction, the number ``Mz`` in
18 the `z` direction within the ice, and the number ``Mbz`` in the `z` direction
19 within the bedrock thermal layer. These are specified by options :opt:`-Mx`, :opt:`-My`,
20 :opt:`-Mz`, and :opt:`-Mbz`, respectively. The defaults are 61, 61, 31, and 1,
21 respectively. Note that ``Mx``, ``My``, ``Mz``, and ``Mbz`` all indicate the number of
22 grid *points* so the number of grid *spaces* are one less, thus 60, 60, 30, and 0 in the
23 default case.
24
25 The lowest grid point in a column of ice, at `z=0`, coincides with the highest grid point
26 in the bedrock, so ``Mbz`` must always be at least one. Choosing ``Mbz`` `>1` is required
27 to use the bedrock thermal model. When a thermal bedrock layer is used, the distance
28 ``Lbz`` is controlled by the :opt:`-Lbz` option. Note that ``Mbz`` is unrelated to the bed
29 deformation model (glacial isostasy model); see section :ref:`sec-beddef`.
30
31 In the quadratically-spaced case the vertical spacing near the ice/bedrock interface is
32 about four times finer than it would be with equal spacing for the same value of ``Mz``,
33 while the spacing near the top of the computational box is correspondingly coarser. For a
34 detailed description of the spacing of the grid, see the documentation on
35 ``IceGrid::compute_vertical_levels()`` in the `PISM class browser <pism-browser_>`_.
36
37 The user should specify the grid when using ``-bootstrap`` or when initializing a
38 verification test (section :ref:`sec-verif`) or a simplified-geometry experiment (section
39 :ref:`sec-simp`). If one initializes PISM from a saved model state using ``-i`` then the
40 input file determines all grid parameters. For instance, the command
41
42 .. code-block:: none
43
44 pismr -i foo.nc -y 100
45
46 should work fine if ``foo.nc`` is a PISM output file. Because ``-i`` input files take
47 precedence over options,
48
49 .. code-block:: none
50
51 pismr -i foo.nc -Mz 201 -y 100
52
53 will give a warning that "``PISM WARNING: ignoring command-line option '-Mz'``".
54
55 .. _sec-grid-registration:
56
57 Grid registration
58 ^^^^^^^^^^^^^^^^^
59
60 .. |xmin| replace:: `x_\text{min}`
61 .. |xmax| replace:: `x_\text{max}`
62
63 PISM's horizontal computational grid is uniform and (by default) cell-centered.\ [#]_
64
65 This is not the only possible interpretation, but it is consistent with the finite-volume
66 handling of mass (ice thickness) evolution is PISM.
67
68 Consider a grid with minimum and maximum `x` coordinates `x_\text{min}` and `x_\text{max}`
69 and the spacing `\Delta x`. The cell-centered interpretation implies that the domain
70 extends *past* |xmin| and |xmax| by one half of the grid spacing, see
71 :numref:`fig-cell-center`.
72
73 .. figure:: figures/grid-centered-both.png
74 :name: fig-cell-center
75 :width: 80%
76 :align: center
77
78 Computational grids using the default (``center``) grid registration.
79
80 *Left*: a coarse grid. *Right*: a finer grid covering the same domain.
81
82 The solid black line represents the domain boundary, dashed red lines are cell
83 boundaries, black circles represent grid points.
84
85 When getting the size of the domain from an input file, PISM will compute grid parameters
86 as follows:
87
88 .. math::
89
90 \Delta x &= x_1 - x_0
91
92 L_x &= \frac12 ((x_\text{max} - x_\text{min}) + \Delta x).
93
94 This is not an issue when re-starting from a PISM output file but can cause confusion when
95 specifying grid parameters at bootstrapping and reading in fields using "regridding."
96
97 For example:
98
99 .. code-block:: bash
100
101 > pisms -grid.registration center \
102 -Lx 10 -Mx 4 \
103 -y 0 -verbose 1 \
104 -o grid-test.nc
105 > ncdump -v x grid-test.nc | tail -2 | head -1
106 x = -7500, -2500, 2500, 7500 ;
107
108 Note that we specified the domain half width of 10 km and selected 4 grid points in the
109 `x` direction. The resulting `x` coordinates range from `-7500` meters to `7500` meters with
110 the grid spacing of `5` km.
111
112 In summary, with the default (center) grid registration
113
114 .. math::
115 :label: eq-grid-center
116
117 \Delta x &= \frac{2 L_x}{M_x},
118
119 x_\text{min} &= x_c - L_x + \frac12 \Delta x,
120
121 x_\text{max} &= x_c + L_x - \frac12 \Delta x,
122
123 where `x_c` is the `x`\-coordinate of the domain center.
124
125 .. note::
126
127 One advantage of this approach is that it is easy to build a set of grids covering a
128 given region such that grid cells nest within each other as in
129 :numref:`fig-cell-center`. In particular, this makes it easier to create a set of
130 surface mass balance fields for the domain that use different resolutions but *have the
131 same total SMB*.
132
133 Compare this to
134
135 .. code-block:: bash
136
137 > pisms -grid.registration corner \
138 -Lx 10 -Mx 5 \
139 -y 0 -verbose 1 \
140 -o grid-test.nc
141 > ncdump -v x grid-test.nc | tail -2 | head -1
142 x = -10000, -5000, 0, 5000, 10000 ;
143
144 Here the grid spacing is also 5 km, although there are 5 grid points in the `x` direction
145 and `x` coordinates range from `-10000` to `10000`.
146
147 With the "corner" grid registration
148
149 .. math::
150 :label: eq-grid-corner
151
152 \Delta x &= \frac{2 L_x}{M_x - 1},
153
154 x_\text{min} &= x_c - L_x,
155
156 x_\text{max} &= x_c + L_x.
157
158 See :numref:`fig-cell-corner` for an illustration.
159
160 .. figure:: figures/grid-corner-both.png
161 :name: fig-cell-corner
162 :width: 80%
163 :align: center
164
165 Computational grids using the ``corner`` grid registration.
166
167 *Left*: a coarse grid. *Right*: a finer grid covering the same domain.
168
169 To switch between :eq:`eq-grid-center` and :eq:`eq-grid-corner`,
170 set the configuration parameter :config:`grid.registration`.
171
172 .. _sec-projections:
173
174 Grid projections
175 ^^^^^^^^^^^^^^^^
176
177 PISM can use the PROJ_ library (see :ref:`sec-install-prerequisites`) and projection
178 information to compute
179
180 - latitudes and longitudes of grid points (variables :var:`lat` and :var:`lon`), and
181 - latitudes and longitudes of cell corners (variables :var:`lat_bnds` and :var:`lon_bnds`).
182
183 To use this feature, compile PISM with PROJ and add the global attribute ``proj``
184 containing the parameter string describing the projection to the input file.
185
186 For example, the input file ``pism_Greenland_5km_v1.1.nc`` in :ref:`sec-start` has the
187 following:
188
189 .. code-block:: bash
190
191 > ncdump -h pism_Greenland_5km_v1.1.nc | grep :proj
192 :proj = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" ;
193
194 The spinup run in that example disables the code re-computing longitude, latitude grid
195 coordinates using projection information to avoid the dependency on PROJ (look for
196 :opt:`-grid.recompute_longitude_and_latitude` in the command). If we remove this option,
197 PISM will report the following.
198
199 .. code-block:: none
200
201 > pismr -i pism_Greenland_5km_v1.1.nc \
202 -bootstrap -Mx 76 -My 141 -Mz 101 -Mbz 11 ... \
203 -grid.recompute_longitude_and_latitude true ... -o output.nc
204 ...
205 * Got projection parameters "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" from "pism_Greenland_5km_v1.1.nc".
206 * Computing longitude and latitude using projection parameters...
207 ...
208 ... done with run
209 Writing model state to file `output.nc'...
210
211 If the ``proj`` attribute contains the string "``+init=epsg:XXXX``" where ``XXXX`` is
212 3413, 3031, or 26710, PISM will also create a CF-conforming ``mapping`` variable
213 describing the projection in use.
214
215 "Mapping" variables following CF metadata conventions in input files are copied to output
216 files (including ``-extra_file``\s) but are **not** used to compute latitude/longitude
217 coordinates.
218
219 To simplify post-processing and analysis with CDO PISM adds the PROJ string (if known)
220 to the mapping variable, putting it in the ``proj_params`` attribute.
221
222 .. code-block:: none
223
224 > ncdump -h g20km_10ka_hy.nc | grep mapping:
225
226 mapping:ellipsoid = "WGS84" ;
227 mapping:grid_mapping_name = "polar_stereographic" ;
228 mapping:false_easting = 0. ;
229 mapping:false_northing = 0. ;
230 mapping:latitude_of_projection_origin = 90. ;
231 mapping:standard_parallel = 71. ;
232 mapping:straight_vertical_longitude_from_pole = -39. ;
233 mapping:proj_params = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" ;
234
235 .. _sec-domain-distribution:
236
237 Parallel domain distribution
238 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
239
240 When running PISM in parallel with ``mpiexec -n N``, the horizontal grid is distributed
241 across `N` processes [#]_. PISM divides the grid into `N_x` parts in the `x` direction and
242 `N_y` parts in the `y` direction. By default this is done automatically, with the goal
243 that `N_x\times N_y = N` and `N_x` is as close to `N_y` as possible. Note that `N` should,
244 therefore, be a composite (not prime) number.
245
246 Users seeking to override this default can specify `N_x` and `N_y` using the :opt:`-Nx`
247 and :opt:`-Ny` command-line options.
248
249 Once `N_x` and `N_y` are computed, PISM computes sizes of sub-domains `M_{x,i}` so that
250 `\sum_{i=1}^{N_x}M_{x,i} = M_x` and `M_{x,i} - \left\lfloor M_x / N_x \right\rfloor < 1`.
251 To specify strip widths `M_{x,i}` and `M_{y,i}`, use command-line options :opt:`-procs_x`
252 and :opt:`-procs_y`. Each option takes a comma-separated list of numbers as its argument.
253 For example,
254
255 .. code-block:: none
256
257 mpiexec -n 3 pisms -Mx 101 -My 101 \
258 -Nx 1 -procs_x 101 \
259 -Ny 3 -procs_y 20,61,20
260
261 splits a `101 \times 101` grid into 3 strips along the `x` axis.
262
263 To see the parallel domain decomposition from a completed run, see the :var:`rank`
264 variable in the output file, e.g. using ``-o_size big``. The same :var:`rank` variable is
265 available as a spatial diagnostic field (section :ref:`sec-saving-diagnostics`).
266
267 .. rubric:: Footnotes
268
269 .. [#] This is consistent with the `CF Conventions`_ document for data-sets without cell
270 bounds: "If bounds are not provided, an application might reasonably assume the
271 gridpoints to be at the centers of the cells, but we do not require that in this
272 standard."
273
274 .. [#] In most cases one process corresponds to one "core" of your computer.