PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
Geometry.cc
Go to the documentation of this file.
1 /* Copyright (C) 2017, 2018, 2019, 2020, 2021 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 "Geometry.hh"
21 
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/IceModelVec2CellType.hh"
24 #include "pism/util/Mask.hh"
25 #include "pism/util/pism_utilities.hh"
26 #include "pism/geometry/grounded_cell_fraction.hh"
27 #include "pism/util/Context.hh"
28 #include "pism/util/VariableMetadata.hh"
29 #include "pism/util/io/File.hh"
30 #include "pism/util/io/io_helpers.hh"
31 
32 namespace pism {
33 
35  // FIXME: ideally these fields should be "global", i.e. without ghosts.
36  // (However this may increase communication costs...)
37  : m_stencil_width(static_cast<int>(grid->ctx()->config()->get_number("grid.max_stencil_width"))),
38  latitude(grid, "lat", WITHOUT_GHOSTS),
39  longitude(grid, "lon", WITHOUT_GHOSTS),
40  bed_elevation(grid, "topg", WITH_GHOSTS, m_stencil_width),
41  sea_level_elevation(grid, "sea_level", WITH_GHOSTS),
42  ice_thickness(grid, "thk", WITH_GHOSTS, m_stencil_width),
43  ice_area_specific_volume(grid, "ice_area_specific_volume", WITH_GHOSTS),
44  cell_type(grid, "mask", WITH_GHOSTS, m_stencil_width),
45  cell_grounded_fraction(grid, "cell_grounded_fraction", WITHOUT_GHOSTS),
46  ice_surface_elevation(grid, "usurf", WITH_GHOSTS, m_stencil_width) {
47 
48  latitude.set_attrs("mapping", "latitude", "degree_north", "degree_north", "latitude", 0);
50  latitude.metadata()["grid_mapping"] = "";
51  latitude.metadata()["valid_range"] = {-90.0, 90.0};
52 
53  longitude.set_attrs("mapping", "longitude", "degree_east", "degree_east", "longitude", 0);
55  longitude.metadata()["grid_mapping"] = "";
56  longitude.metadata()["valid_range"] = {-180.0, 180.0};
57 
58  bed_elevation.set_attrs("model_state", "bedrock surface elevation",
59  "m", "m", "bedrock_altitude", 0);
60 
61  sea_level_elevation.set_attrs("model_state",
62  "sea level elevation above reference ellipsoid", "meters", "meters",
63  "sea_surface_height_above_reference_ellipsoid", 0);
64 
65  ice_thickness.set_attrs("model_state", "land ice thickness",
66  "m", "m", "land_ice_thickness", 0);
67  ice_thickness.metadata()["valid_min"] = {0.0};
68 
69  ice_area_specific_volume.set_attrs("model_state",
70  "ice-volume-per-area in partially-filled grid cells",
71  "m3/m2", "m3/m2", "", 0);
72  ice_area_specific_volume.metadata()["comment"] =
73  "this variable represents the amount of ice in a partially-filled cell and not "
74  "the corresponding geometry, so thinking about it as 'thickness' is not helpful";
75 
76  cell_type.set_attrs("diagnostic", "ice-type (ice-free/grounded/floating/ocean) integer mask",
77  "", "", "", 0);
79  cell_type.metadata()["flag_meanings"] =
80  "ice_free_bedrock grounded_ice floating_ice ice_free_ocean";
82 
84  "fractional grounded/floating mask (floating=0, grounded=1)",
85  "", "", "", 0);
86 
87  ice_surface_elevation.set_attrs("diagnostic", "ice upper surface elevation",
88  "m", "m", "surface_altitude", 0);
89 
90  // make sure all the fields are initialized
91  latitude.set(0.0);
92  longitude.set(0.0);
93  bed_elevation.set(0.0);
95  ice_thickness.set(0.0);
97  ensure_consistency(0.0);
98 }
99 
100 void check_minimum_ice_thickness(const IceModelVec2S &ice_thickness) {
101  IceGrid::ConstPtr grid = ice_thickness.grid();
102 
103  IceModelVec::AccessList list(ice_thickness);
104 
105  ParallelSection loop(grid->com);
106  try {
107  for (Points p(*grid); p; p.next()) {
108  const int i = p.i(), j = p.j();
109 
110  if (ice_thickness(i, j) < 0.0) {
112  "H = %e (negative) at point i=%d, j=%d",
113  ice_thickness(i, j), i, j);
114  }
115  }
116  } catch (...) {
117  loop.failed();
118  }
119  loop.check();
120 }
121 
122 void Geometry::ensure_consistency(double ice_free_thickness_threshold) {
124  Config::ConstPtr config = grid->ctx()->config();
125 
127 
131 
132  // first ensure that ice_area_specific_volume is 0 if ice_thickness > 0.
133  {
134  ParallelSection loop(grid->com);
135  try {
136  for (Points p(*grid); p; p.next()) {
137  const int i = p.i(), j = p.j();
138 
139  if (ice_thickness(i, j) > 0.0 and ice_area_specific_volume(i, j) > 0.0) {
141  ice_area_specific_volume(i, j) = 0.0;
142  }
143  }
144  } catch (...) {
145  loop.failed();
146  }
147  loop.check();
148  }
149 
150  // compute cell type and surface elevation
151  {
152  GeometryCalculator gc(*config);
153  gc.set_icefree_thickness(ice_free_thickness_threshold);
154 
155  ParallelSection loop(grid->com);
156  try {
157  for (Points p(*grid); p; p.next()) {
158  const int i = p.i(), j = p.j();
159 
160  int mask = 0;
162  &mask, &ice_surface_elevation(i, j));
163  cell_type(i, j) = mask;
164  }
165  } catch (...) {
166  loop.failed();
167  }
168  loop.check();
169  }
170 
175 
176  const double
177  ice_density = config->get_number("constants.ice.density"),
178  ocean_density = config->get_number("constants.sea_water.density");
179 
180  try {
181  compute_grounded_cell_fraction(ice_density,
182  ocean_density,
187  } catch (RuntimeError &e) {
188  e.add_context("computing the grounded cell fraction");
189 
190  std::string output_file = config->get_string("output.file_name");
191  std::string o_file = filename_add_suffix(output_file,
192  "_grounded_cell_fraction_failed", "");
193  // save geometry to a file for debugging
194  dump(o_file.c_str());
195  throw;
196  }
197 }
198 
199 void Geometry::dump(const char *filename) const {
200  auto grid = ice_thickness.grid();
201 
202  File file(grid->com, filename,
203  string_to_backend(grid->ctx()->config()->get_string("output.format")),
205  grid->ctx()->pio_iosys_id());
206 
207  io::define_time(file, *grid->ctx());
208  io::append_time(file, *grid->ctx()->config(), 0.0);
209 
210  latitude.write(file);
211  longitude.write(file);
212  bed_elevation.write(file);
214  ice_thickness.write(file);
216  cell_type.write(file);
219 }
220 
221 /*! Compute the elevation of the bottom surface of the ice.
222  */
223 void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result) {
224 
225  auto grid = result.grid();
226  auto config = grid->ctx()->config();
227 
228  double
229  ice_density = config->get_number("constants.ice.density"),
230  water_density = config->get_number("constants.sea_water.density"),
231  alpha = ice_density / water_density;
232 
233  const IceModelVec2S &ice_thickness = geometry.ice_thickness;
234  const IceModelVec2S &bed_elevation = geometry.bed_elevation;
235  const IceModelVec2S &sea_level = geometry.sea_level_elevation;
236 
237  IceModelVec::AccessList list{&ice_thickness, &bed_elevation, &sea_level, &result};
238 
239  ParallelSection loop(grid->com);
240  try {
241  for (Points p(*grid); p; p.next()) {
242  const int i = p.i(), j = p.j();
243 
244  double
245  b_grounded = bed_elevation(i, j),
246  b_floating = sea_level(i, j) - alpha * ice_thickness(i, j);
247 
248  result(i, j) = std::max(b_grounded, b_floating);
249  }
250  } catch (...) {
251  loop.failed();
252  }
253  loop.check();
254 
255  result.update_ghosts();
256 }
257 
258 //! Computes the ice volume, in m^3.
259 double ice_volume(const Geometry &geometry, double thickness_threshold) {
260  auto grid = geometry.ice_thickness.grid();
261  auto config = grid->ctx()->config();
262 
263  IceModelVec::AccessList list{&geometry.ice_thickness};
264 
265  double volume = 0.0;
266 
267  auto cell_area = grid->cell_area();
268 
269  {
270  for (Points p(*grid); p; p.next()) {
271  const int i = p.i(), j = p.j();
272 
273  if (geometry.ice_thickness(i,j) >= thickness_threshold) {
274  volume += geometry.ice_thickness(i,j) * cell_area;
275  }
276  }
277  }
278 
279  // Add the volume of the ice in Href:
280  if (config->get_flag("geometry.part_grid.enabled")) {
281  list.add(geometry.ice_area_specific_volume);
282  for (Points p(*grid); p; p.next()) {
283  const int i = p.i(), j = p.j();
284 
285  volume += geometry.ice_area_specific_volume(i,j) * cell_area;
286  }
287  }
288 
289  return GlobalSum(grid->com, volume);
290 }
291 
293  double thickness_threshold) {
294  auto grid = geometry.ice_thickness.grid();
295  auto config = grid->ctx()->config();
296 
297  const double
298  sea_water_density = config->get_number("constants.sea_water.density"),
299  ice_density = config->get_number("constants.ice.density"),
300  cell_area = grid->cell_area();
301 
302  IceModelVec::AccessList list{&geometry.cell_type, &geometry.ice_thickness,
303  &geometry.bed_elevation, &geometry.sea_level_elevation};
304 
305  double volume = 0.0;
306 
307  for (Points p(*grid); p; p.next()) {
308  const int i = p.i(), j = p.j();
309 
310  const double
311  bed = geometry.bed_elevation(i, j),
312  thickness = geometry.ice_thickness(i, j),
313  sea_level = geometry.sea_level_elevation(i, j);
314 
315  if (geometry.cell_type.grounded(i, j) and thickness > thickness_threshold) {
316  const double cell_ice_volume = thickness * cell_area;
317  if (bed > sea_level) {
318  volume += cell_ice_volume;
319  } else {
320  const double max_floating_volume = (sea_level - bed) * cell_area * (sea_water_density / ice_density);
321  volume += cell_ice_volume - max_floating_volume;
322  }
323  }
324  } // end of the loop over grid points
325 
326  return GlobalSum(grid->com, volume);
327 }
328 
329 //! Computes ice area, in m^2.
330 double ice_area(const Geometry &geometry, double thickness_threshold) {
331  auto grid = geometry.ice_thickness.grid();
332 
333  double area = 0.0;
334 
335  auto cell_area = grid->cell_area();
336 
337  IceModelVec::AccessList list{&geometry.ice_thickness};
338  for (Points p(*grid); p; p.next()) {
339  const int i = p.i(), j = p.j();
340 
341  if (geometry.ice_thickness(i, j) >= thickness_threshold) {
342  area += cell_area;
343  }
344  }
345 
346  return GlobalSum(grid->com, area);
347 }
348 
349 //! Computes grounded ice area, in m^2.
350 double ice_area_grounded(const Geometry &geometry, double thickness_threshold) {
351  auto grid = geometry.ice_thickness.grid();
352 
353  double area = 0.0;
354 
355  auto cell_area = grid->cell_area();
356 
357  IceModelVec::AccessList list{&geometry.cell_type, &geometry.ice_thickness};
358  for (Points p(*grid); p; p.next()) {
359  const int i = p.i(), j = p.j();
360 
361  if (geometry.cell_type.grounded(i, j) and
362  geometry.ice_thickness(i, j) >= thickness_threshold) {
363  area += cell_area;
364  }
365  }
366 
367  return GlobalSum(grid->com, area);
368 }
369 
370 //! Computes floating ice area, in m^2.
371 double ice_area_floating(const Geometry &geometry, double thickness_threshold) {
372  auto grid = geometry.ice_thickness.grid();
373 
374  double area = 0.0;
375 
376  auto cell_area = grid->cell_area();
377 
378  IceModelVec::AccessList list{&geometry.cell_type, &geometry.ice_thickness};
379  for (Points p(*grid); p; p.next()) {
380  const int i = p.i(), j = p.j();
381 
382  if (geometry.cell_type.ocean(i, j) and
383  geometry.ice_thickness(i, j) >= thickness_threshold) {
384  area += cell_area;
385  }
386  }
387 
388  return GlobalSum(grid->com, area);
389 }
390 
391 
392 //! Computes the sea level rise that would result if all the ice were melted.
393 double sea_level_rise_potential(const Geometry &geometry, double thickness_threshold) {
394  auto config = geometry.ice_thickness.grid()->ctx()->config();
395 
396  const double
397  water_density = config->get_number("constants.fresh_water.density"),
398  ice_density = config->get_number("constants.ice.density"),
399  ocean_area = config->get_number("constants.global_ocean_area");
400 
401  const double
402  volume = ice_volume_not_displacing_seawater(geometry,
403  thickness_threshold),
404  additional_water_volume = (ice_density / water_density) * volume,
405  sea_level_change = additional_water_volume / ocean_area;
406 
407  return sea_level_change;
408 }
409 
410 
411 /*!
412  * @brief Set no_model_mask variable to have value 1 in strip of width 'strip' m around
413  * edge of computational domain, and value 0 otherwise.
414  */
415 void set_no_model_strip(const IceGrid &grid, double width, IceModelVec2Int &result) {
416 
417  if (width <= 0.0) {
418  return;
419  }
420 
421  IceModelVec::AccessList list(result);
422 
423  for (Points p(grid); p; p.next()) {
424  const int i = p.i(), j = p.j();
425 
426  if (in_null_strip(grid, i, j, width)) {
427  result(i, j) = 1;
428  } else {
429  result(i, j) = 0;
430  }
431  }
432 
433  result.update_ghosts();
434 }
435 
436 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
std::shared_ptr< const Config > ConstPtr
High-level PISM I/O class.
Definition: File.hh:51
void compute(const IceModelVec2S &sea_level, const IceModelVec2S &bed, const IceModelVec2S &thickness, IceModelVec2Int &out_mask, IceModelVec2S &out_surface) const
Definition: Mask.cc:26
void set_icefree_thickness(double threshold)
Definition: Mask.hh:79
Geometry(const IceGrid::ConstPtr &grid)
Definition: Geometry.cc:34
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
IceModelVec2S latitude
Definition: Geometry.hh:43
IceModelVec2S longitude
Definition: Geometry.hh:44
void ensure_consistency(double ice_free_thickness_threshold)
Definition: Geometry.cc:122
IceModelVec2S cell_grounded_fraction
Definition: Geometry.hh:58
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S ice_area_specific_volume
Definition: Geometry.hh:54
void dump(const char *filename) const
Definition: Geometry.cc:199
IceModelVec2S sea_level_elevation
Definition: Geometry.hh:50
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
Describes the PISM grid and the distribution of data across processors.
Definition: IceGrid.hh:228
bool ocean(int i, int j) const
bool grounded(int i, int j) const
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
void add(double alpha, const IceModelVec2S &x)
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
void set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
Definition: iceModelVec.cc:175
void failed()
Indicates a failure of a parallel section.
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
void set_output_type(IO_Type type)
#define PISM_ERROR_LOCATION
void append_time(const File &file, const Config &config, double time_seconds)
Prepare a file for output.
Definition: io_helpers.cc:257
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
Definition: io_helpers.cc:220
double ice_volume_not_displacing_seawater(const Geometry &geometry, double thickness_threshold)
Definition: Geometry.cc:292
double ice_area(const Geometry &geometry, double thickness_threshold)
Computes ice area, in m^2.
Definition: Geometry.cc:330
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
@ PISM_INT
Definition: IO_Flags.hh:36
void compute_grounded_cell_fraction(double ice_density, double ocean_density, const IceModelVec2S &sea_level, const IceModelVec2S &ice_thickness, const IceModelVec2S &bed_topography, IceModelVec2S &result)
IO_Backend string_to_backend(const std::string &backend)
Definition: File.cc:60
double ice_area_floating(const Geometry &geometry, double thickness_threshold)
Computes floating ice area, in m^2.
Definition: Geometry.cc:371
double sea_level_rise_potential(const Geometry &geometry, double thickness_threshold)
Computes the sea level rise that would result if all the ice were melted.
Definition: Geometry.cc:393
@ MASK_FLOATING
Definition: Mask.hh:33
@ MASK_ICE_FREE_OCEAN
Definition: Mask.hh:34
@ MASK_ICE_FREE_BEDROCK
Definition: Mask.hh:31
@ MASK_GROUNDED
Definition: Mask.hh:32
std::string filename_add_suffix(const std::string &filename, const std::string &separator, const std::string &suffix)
Adds a suffix to a filename.
double ice_volume(const Geometry &geometry, double thickness_threshold)
Computes the ice volume, in m^3.
Definition: Geometry.cc:259
void set_no_model_strip(const IceGrid &grid, double width, IceModelVec2Int &result)
Set no_model_mask variable to have value 1 in strip of width 'strip' m around edge of computational d...
Definition: Geometry.cc:415
@ PISM_READWRITE_CLOBBER
create a file for writing, overwrite if present
Definition: IO_Flags.hh:53
void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result)
Definition: Geometry.cc:223
double ice_area_grounded(const Geometry &geometry, double thickness_threshold)
Computes grounded ice area, in m^2.
Definition: Geometry.cc:350
bool in_null_strip(const IceGrid &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: IceGrid.hh:339
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
void check_minimum_ice_thickness(const IceModelVec2S &ice_thickness)
Definition: Geometry.cc:100
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49