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