tFrontRetreat.cc - 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
---
tFrontRetreat.cc (9632B)
---
1 /* Copyright (C) 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 #include "FrontRetreat.hh"
20
21 #include "pism/util/iceModelVec.hh"
22 #include "pism/util/IceModelVec2CellType.hh"
23 #include "pism/util/MaxTimestep.hh"
24 #include "pism/util/pism_utilities.hh"
25 #include "pism/geometry/part_grid_threshold_thickness.hh"
26 #include "pism/geometry/Geometry.hh"
27
28 namespace pism {
29
30 FrontRetreat::FrontRetreat(IceGrid::ConstPtr g)
31 : Component(g),
32 m_tmp(m_grid, "temporary_storage", WITH_GHOSTS, 1) {
33
34 m_tmp.set_attrs("internal", "additional mass loss at points near the front",
35 "m", "m", "", 0);
36
37 m_cell_type.create(m_grid, "cell_type", WITH_GHOSTS, 1);
38 m_cell_type.set_attrs("internal", "cell type mask", "", "", "", 0);
39 }
40
41 FrontRetreat::~FrontRetreat() {
42 // empty
43 }
44
45 /*!
46 * Compute the modified mask to avoid "wrapping around" of front retreat at domain
47 * boundaries.
48 */
49 void FrontRetreat::compute_modified_mask(const IceModelVec2CellType &input,
50 IceModelVec2CellType &output) const {
51
52 IceModelVec::AccessList list{&input, &output};
53
54 const int Mx = m_grid->Mx();
55 const int My = m_grid->My();
56
57 ParallelSection loop(m_grid->com);
58 try {
59 for (PointsWithGhosts p(*m_grid); p; p.next()) {
60 const int i = p.i(), j = p.j();
61
62 if (i < 0 or i >= Mx or j < 0 or j >= My) {
63 output(i, j) = MASK_ICE_FREE_OCEAN;
64 } else {
65 output(i, j) = input(i, j);
66 }
67 }
68 } catch (...) {
69 loop.failed();
70 }
71 loop.check();
72 }
73
74 /*!
75 * Compute the maximum time step length provided a horizontal retreat rate.
76 */
77 MaxTimestep FrontRetreat::max_timestep(const IceModelVec2CellType &cell_type,
78 const IceModelVec2Int &bc_mask,
79 const IceModelVec2S &retreat_rate) const {
80
81 IceGrid::ConstPtr grid = retreat_rate.grid();
82 units::System::Ptr sys = grid->ctx()->unit_system();
83
84 using units::convert;
85
86 // About 9 hours which corresponds to 10000 km year-1 on a 10 km grid
87 double dt_min = convert(sys, 0.001, "years", "seconds");
88
89 double
90 retreat_rate_max = 0.0,
91 retreat_rate_mean = 0.0;
92 int N_cells = 0;
93
94 IceModelVec::AccessList list{&cell_type, &bc_mask, &retreat_rate};
95
96 for (Points pt(*grid); pt; pt.next()) {
97 const int i = pt.i(), j = pt.j();
98
99 if (cell_type.ice_free_ocean(i, j) and
100 cell_type.next_to_ice(i, j) and
101 bc_mask(i, j) < 0.5) {
102 // NB: this condition has to match the one in update_geometry()
103
104 const double C = retreat_rate(i, j);
105
106 N_cells += 1;
107 retreat_rate_mean += C;
108 retreat_rate_max = std::max(C, retreat_rate_max);
109 }
110 }
111
112 N_cells = GlobalSum(grid->com, N_cells);
113 retreat_rate_mean = GlobalSum(grid->com, retreat_rate_mean);
114 retreat_rate_max = GlobalMax(grid->com, retreat_rate_max);
115
116 if (N_cells > 0.0) {
117 retreat_rate_mean /= N_cells;
118 } else {
119 retreat_rate_mean = 0.0;
120 }
121
122 double denom = retreat_rate_max / grid->dx();
123 const double epsilon = convert(sys, 0.001 / (grid->dx() + grid->dy()), "seconds", "years");
124
125 double dt = 1.0 / (denom + epsilon);
126
127 m_log->message(3,
128 " frontal retreat: maximum rate = %.2f m/year gives dt=%.5f years\n"
129 " mean rate = %.2f m/year over %d cells\n",
130 convert(m_sys, retreat_rate_max, "m second-1", "m year-1"),
131 convert(m_sys, dt, "seconds", "years"),
132 convert(m_sys, retreat_rate_mean, "m second-1", "m year-1"),
133 N_cells);
134
135 return MaxTimestep(std::max(dt, dt_min), "front_retreat");
136 }
137
138 /*!
139 * Update ice geometry by applying a horizontal retreat rate.
140 *
141 * This code applies a horizontal retreat rate at "partially-filled" cells that are next
142 * to icy cells.
143 *
144 * Models providing the "retreat rate" should set this field to zero in areas where a
145 * particular parameterization does not apply. (For example: some calving models apply at
146 * shelf calving fronts, others may apply at grounded termini but not at ice shelves,
147 * etc).
148 */
149 void FrontRetreat::update_geometry(double dt,
150 const Geometry &geometry,
151 const IceModelVec2Int &bc_mask,
152 const IceModelVec2S &retreat_rate,
153 IceModelVec2S &Href,
154 IceModelVec2S &ice_thickness) {
155
156 const IceModelVec2S &bed = geometry.bed_elevation;
157 const IceModelVec2S &sea_level = geometry.sea_level_elevation;
158 const IceModelVec2S &surface_elevation = geometry.ice_surface_elevation;
159
160 if (m_config->get_flag("geometry.front_retreat.wrap_around")) {
161 m_cell_type.copy_from(geometry.cell_type);
162 } else {
163 // compute the modified mask needed to prevent front retreat from "wrapping around"
164 // the computational domain
165 compute_modified_mask(geometry.cell_type, m_cell_type);
166 }
167
168 const double dx = m_grid->dx();
169
170 m_tmp.set(0.0);
171
172 IceModelVec::AccessList list{&ice_thickness, &bc_mask,
173 &bed, &sea_level, &m_cell_type, &Href, &m_tmp, &retreat_rate,
174 &surface_elevation};
175
176 // Prepare to loop over neighbors: directions
177 const Direction dirs[] = {North, East, South, West};
178
179 // Step 1: Apply the computed horizontal retreat rate:
180 for (Points pt(*m_grid); pt; pt.next()) {
181 const int i = pt.i(), j = pt.j();
182
183 // apply retreat rate at the margin (i.e. to partially-filled cells) only
184 if (m_cell_type.ice_free_ocean(i, j) and
185 m_cell_type.next_to_ice(i, j) and
186 bc_mask(i, j) < 0.5) {
187 // NB: this condition has to match the one in max_timestep()
188
189 const double
190 rate = retreat_rate(i, j),
191 Href_old = Href(i, j);
192
193 // Compute the number of floating neighbors and the neighbor-averaged ice thickness:
194 double H_threshold = part_grid_threshold_thickness(m_cell_type.int_star(i, j),
195 ice_thickness.star(i, j),
196 surface_elevation.star(i, j),
197 bed(i, j));
198
199 // Calculate mass loss with respect to the associated ice thickness and the grid size:
200 const double Href_change = -dt * rate * H_threshold / dx; // in m
201
202 if (Href_old + Href_change >= 0.0) {
203 // Href is high enough to absorb the mass loss
204 Href(i, j) = Href_old + Href_change;
205 } else {
206 Href(i, j) = 0.0;
207 // Href + Href_change is negative: need to distribute mass loss to neighboring points
208
209 // Find the number of neighbors to distribute to.
210 //
211 // We consider floating cells and grounded cells with the base below sea level. In other
212 // words, additional mass losses are distributed to shelf calving fronts and grounded marine
213 // termini.
214 int N = 0;
215 {
216 auto
217 M = m_cell_type.int_star(i, j),
218 BC = bc_mask.int_star(i, j);
219
220 auto
221 b = bed.star(i, j),
222 sl = sea_level.star(i, j);
223
224 for (int n = 0; n < 4; ++n) {
225 Direction direction = dirs[n];
226 int m = M[direction];
227 int bc = BC[direction];
228
229 // note: this is where the modified cell type mask is used to prevent
230 // wrap-around
231 if (bc == 0 and // distribute to regular (*not* Dirichlet B.C.) neighbors only
232 (mask::floating_ice(m) or
233 (mask::grounded_ice(m) and b[direction] < sl[direction]))) {
234 N += 1;
235 }
236 }
237 }
238
239 if (N > 0) {
240 m_tmp(i, j) = (Href_old + Href_change) / (double)N;
241 } else {
242 // No shelf calving front of grounded terminus to distribute to: retreat stops here.
243 m_tmp(i, j) = 0.0;
244 }
245 }
246
247 } // end of "if ice free ocean next to ice and not a BC location "
248 } // end of loop over grid points
249
250 // Step 2: update ice thickness and Href in neighboring cells if we need to propagate mass losses.
251 m_tmp.update_ghosts();
252
253 for (Points p(*m_grid); p; p.next()) {
254 const int i = p.i(), j = p.j();
255
256 // Note: this condition has to match the one in step 1 above.
257 if (bc_mask.as_int(i, j) == 0 and
258 (m_cell_type.floating_ice(i, j) or
259 (m_cell_type.grounded_ice(i, j) and bed(i, j) < sea_level(i, j)))) {
260
261 const double delta_H = (m_tmp(i + 1, j) + m_tmp(i - 1, j) +
262 m_tmp(i, j + 1) + m_tmp(i, j - 1));
263
264 if (delta_H < 0.0) {
265 Href(i, j) = ice_thickness(i, j) + delta_H; // in m
266 ice_thickness(i, j) = 0.0;
267 }
268
269 // Stop retreat if the current cell does not have enough ice to absorb the loss.
270 if (Href(i, j) < 0.0) {
271 Href(i, j) = 0.0;
272 }
273 }
274 }
275 }
276
277 } // end of namespace pism