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