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