ShallowStressBalance.cc (12804B)
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Constantine Khroulev and Ed Bueler 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 "ShallowStressBalance.hh" 20 #include "pism/basalstrength/basal_resistance.hh" 21 #include "pism/rheology/FlowLawFactory.hh" 22 23 #include "pism/util/Mask.hh" 24 #include "pism/util/Vars.hh" 25 #include "pism/util/error_handling.hh" 26 #include "pism/util/pism_options.hh" 27 #include "pism/util/IceModelVec2CellType.hh" 28 29 #include "SSB_diagnostics.hh" 30 31 namespace pism { 32 namespace stressbalance { 33 34 //! Evaluate the margin pressure difference term in the calving-front BC. 35 // 36 // Units: (kg / m3) * (m / s2) * m2 = Pa m 37 double margin_pressure_difference(bool shelf, bool dry_mode, double H, double bed, 38 double sea_level, double rho_ice, double rho_ocean, 39 double g) { 40 if (shelf) { 41 // floating shelf 42 return 0.5 * rho_ice * g * (1.0 - (rho_ice / rho_ocean)) * H * H; 43 } else { 44 // grounded terminus 45 if (bed >= sea_level or dry_mode) { 46 return 0.5 * rho_ice * g * H * H; 47 } else { 48 return 0.5 * rho_ice * g * (H * H - (rho_ocean / rho_ice) * pow(sea_level - bed, 2.0)); 49 } 50 } 51 } 52 53 using pism::mask::ice_free; 54 55 ShallowStressBalance::ShallowStressBalance(IceGrid::ConstPtr g) 56 : Component(g), m_basal_sliding_law(NULL), m_flow_law(NULL), m_EC(g->ctx()->enthalpy_converter()) { 57 58 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width"); 59 60 if (m_config->get_flag("basal_resistance.pseudo_plastic.enabled") == true) { 61 m_basal_sliding_law = new IceBasalResistancePseudoPlasticLaw(*m_config); 62 } else { 63 m_basal_sliding_law = new IceBasalResistancePlasticLaw(*m_config); 64 } 65 66 m_velocity.create(m_grid, "bar", WITH_GHOSTS, WIDE_STENCIL); // components ubar, vbar 67 m_velocity.set_attrs("model_state", 68 "thickness-advective ice velocity (x-component)", 69 "m s-1", "m s-1", "", 0); 70 m_velocity.set_attrs("model_state", 71 "thickness-advective ice velocity (y-component)", 72 "m s-1", "m s-1", "", 1); 73 74 m_basal_frictional_heating.create(m_grid, "bfrict", WITHOUT_GHOSTS); 75 m_basal_frictional_heating.set_attrs("diagnostic", 76 "basal frictional heating", 77 "W m-2", "mW m-2", "", 0); 78 } 79 80 ShallowStressBalance::~ShallowStressBalance() { 81 delete m_basal_sliding_law; 82 } 83 84 void ShallowStressBalance::init() { 85 this->init_impl(); 86 } 87 88 void ShallowStressBalance::init_impl() { 89 // empty 90 } 91 92 std::string ShallowStressBalance::stdout_report() const { 93 return ""; 94 } 95 96 std::shared_ptr<const rheology::FlowLaw> ShallowStressBalance::flow_law() const { 97 return m_flow_law; 98 } 99 100 EnthalpyConverter::Ptr ShallowStressBalance::enthalpy_converter() const { 101 return m_EC; 102 } 103 104 const IceBasalResistancePlasticLaw* ShallowStressBalance::sliding_law() const { 105 return m_basal_sliding_law; 106 } 107 108 //! \brief Get the thickness-advective 2D velocity. 109 const IceModelVec2V& ShallowStressBalance::velocity() const { 110 return m_velocity; 111 } 112 113 //! \brief Get the basal frictional heating (for the adaptive energy time-stepping). 114 const IceModelVec2S& ShallowStressBalance::basal_frictional_heating() { 115 return m_basal_frictional_heating; 116 } 117 118 119 DiagnosticList ShallowStressBalance::diagnostics_impl() const { 120 DiagnosticList result = { 121 {"beta", Diagnostic::Ptr(new SSB_beta(this))}, 122 {"taub", Diagnostic::Ptr(new SSB_taub(this))}, 123 {"taub_mag", Diagnostic::Ptr(new SSB_taub_mag(this))}, 124 {"taud", Diagnostic::Ptr(new SSB_taud(this))}, 125 {"taud_mag", Diagnostic::Ptr(new SSB_taud_mag(this))} 126 }; 127 128 if(m_config->get_flag("output.ISMIP6")) { 129 result["strbasemag"] = Diagnostic::Ptr(new SSB_taub_mag(this)); 130 } 131 132 return result; 133 } 134 135 136 ZeroSliding::ZeroSliding(IceGrid::ConstPtr g) 137 : ShallowStressBalance(g) { 138 139 // Use the SIA flow law. 140 rheology::FlowLawFactory ice_factory("stress_balance.sia.", m_config, m_EC); 141 m_flow_law = ice_factory.create(); 142 } 143 144 ZeroSliding::~ZeroSliding() { 145 // empty 146 } 147 148 //! \brief Update the trivial shallow stress balance object. 149 void ZeroSliding::update(const Inputs &inputs, bool full_update) { 150 (void) inputs; 151 152 if (full_update) { 153 m_velocity.set(0.0); 154 m_basal_frictional_heating.set(0.0); 155 } 156 } 157 158 //! \brief Compute the basal frictional heating. 159 /*! 160 Ice shelves have zero basal friction heating. 161 162 \param[in] V *basal* sliding velocity 163 \param[in] tauc basal yield stress 164 \param[in] mask (used to determine if floating or grounded) 165 \param[out] result 166 */ 167 void ShallowStressBalance::compute_basal_frictional_heating(const IceModelVec2V &V, 168 const IceModelVec2S &tauc, 169 const IceModelVec2CellType &mask, 170 IceModelVec2S &result) const { 171 172 IceModelVec::AccessList list{&V, &result, &tauc, &mask}; 173 174 for (Points p(*m_grid); p; p.next()) { 175 const int i = p.i(), j = p.j(); 176 177 if (mask.ocean(i,j)) { 178 result(i,j) = 0.0; 179 } else { 180 const double 181 C = m_basal_sliding_law->drag(tauc(i,j), V(i,j).u, V(i,j).v), 182 basal_stress_x = - C * V(i,j).u, 183 basal_stress_y = - C * V(i,j).v; 184 result(i,j) = - basal_stress_x * V(i,j).u - basal_stress_y * V(i,j).v; 185 } 186 } 187 } 188 189 190 SSB_taud::SSB_taud(const ShallowStressBalance *m) 191 : Diag<ShallowStressBalance>(m) { 192 193 // set metadata: 194 m_vars = {SpatialVariableMetadata(m_sys, "taud_x"), 195 SpatialVariableMetadata(m_sys, "taud_y")}; 196 197 set_attrs("X-component of the driving shear stress at the base of ice", "", 198 "Pa", "Pa", 0); 199 set_attrs("Y-component of the driving shear stress at the base of ice", "", 200 "Pa", "Pa", 1); 201 202 for (auto &v : m_vars) { 203 v.set_string("comment", 204 "this field is purely diagnostic (not used by the model)"); 205 } 206 } 207 208 /*! 209 * The driving stress computed here is not used by the model, so this 210 * implementation intentionally does not use the eta-transformation or special 211 * cases at ice margins. 212 */ 213 IceModelVec::Ptr SSB_taud::compute_impl() const { 214 215 IceModelVec2V::Ptr result(new IceModelVec2V); 216 result->create(m_grid, "result", WITHOUT_GHOSTS); 217 result->metadata(0) = m_vars[0]; 218 result->metadata(1) = m_vars[1]; 219 220 const IceModelVec2S *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness"); 221 const IceModelVec2S *surface = m_grid->variables().get_2d_scalar("surface_altitude"); 222 223 double standard_gravity = m_config->get_number("constants.standard_gravity"), 224 ice_density = m_config->get_number("constants.ice.density"); 225 226 IceModelVec::AccessList list{surface, thickness, result.get()}; 227 228 for (Points p(*m_grid); p; p.next()) { 229 const int i = p.i(), j = p.j(); 230 231 double pressure = ice_density * standard_gravity * (*thickness)(i,j); 232 if (pressure <= 0.0) { 233 (*result)(i,j).u = 0.0; 234 (*result)(i,j).v = 0.0; 235 } else { 236 (*result)(i,j).u = - pressure * surface->diff_x_p(i,j); 237 (*result)(i,j).v = - pressure * surface->diff_y_p(i,j); 238 } 239 } 240 241 return result; 242 } 243 244 SSB_taud_mag::SSB_taud_mag(const ShallowStressBalance *m) 245 : Diag<ShallowStressBalance>(m) { 246 247 // set metadata: 248 m_vars = {SpatialVariableMetadata(m_sys, "taud_mag")}; 249 250 set_attrs("magnitude of the gravitational driving stress at the base of ice", "", 251 "Pa", "Pa", 0); 252 m_vars[0].set_string("comment", 253 "this field is purely diagnostic (not used by the model)"); 254 } 255 256 IceModelVec::Ptr SSB_taud_mag::compute_impl() const { 257 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "taud_mag", WITHOUT_GHOSTS)); 258 result->metadata(0) = m_vars[0]; 259 260 IceModelVec2V::Ptr taud = IceModelVec2V::ToVector(SSB_taud(model).compute()); 261 262 result->set_to_magnitude(*taud); 263 264 return result; 265 } 266 267 SSB_taub::SSB_taub(const ShallowStressBalance *m) 268 : Diag<ShallowStressBalance>(m) { 269 // set metadata: 270 m_vars = {SpatialVariableMetadata(m_sys, "taub_x"), 271 SpatialVariableMetadata(m_sys, "taub_y")}; 272 273 set_attrs("X-component of the shear stress at the base of ice", "", 274 "Pa", "Pa", 0); 275 set_attrs("Y-component of the shear stress at the base of ice", "", 276 "Pa", "Pa", 1); 277 278 for (auto &v : m_vars) { 279 v.set_string("comment", 280 "this field is purely diagnostic (not used by the model)"); 281 } 282 } 283 284 285 IceModelVec::Ptr SSB_taub::compute_impl() const { 286 287 IceModelVec2V::Ptr result(new IceModelVec2V); 288 result->create(m_grid, "result", WITHOUT_GHOSTS); 289 result->metadata() = m_vars[0]; 290 result->metadata(1) = m_vars[1]; 291 292 const IceModelVec2V &velocity = model->velocity(); 293 const IceModelVec2S *tauc = m_grid->variables().get_2d_scalar("tauc"); 294 const IceModelVec2CellType &mask = *m_grid->variables().get_2d_cell_type("mask"); 295 296 const IceBasalResistancePlasticLaw *basal_sliding_law = model->sliding_law(); 297 298 IceModelVec::AccessList list{tauc, &velocity, &mask, result.get()}; 299 for (Points p(*m_grid); p; p.next()) { 300 const int i = p.i(), j = p.j(); 301 302 if (mask.grounded_ice(i,j)) { 303 double beta = basal_sliding_law->drag((*tauc)(i,j), velocity(i,j).u, velocity(i,j).v); 304 (*result)(i,j).u = - beta * velocity(i,j).u; 305 (*result)(i,j).v = - beta * velocity(i,j).v; 306 } else { 307 (*result)(i,j).u = 0.0; 308 (*result)(i,j).v = 0.0; 309 } 310 } 311 312 return result; 313 } 314 315 SSB_taub_mag::SSB_taub_mag(const ShallowStressBalance *m) 316 : Diag<ShallowStressBalance>(m) { 317 318 auto ismip6 = m_config->get_flag("output.ISMIP6"); 319 320 // set metadata: 321 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "strbasemag" : "taub_mag")}; 322 323 set_attrs("magnitude of the basal shear stress at the base of ice", 324 "land_ice_basal_drag", // ISMIP6 "standard" name 325 "Pa", "Pa", 0); 326 m_vars[0].set_string("comment", 327 "this field is purely diagnostic (not used by the model)"); 328 } 329 330 IceModelVec::Ptr SSB_taub_mag::compute_impl() const { 331 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "taub_mag", WITHOUT_GHOSTS)); 332 result->metadata(0) = m_vars[0]; 333 334 IceModelVec2V::Ptr taub = IceModelVec2V::ToVector(SSB_taub(model).compute()); 335 336 result->set_to_magnitude(*taub); 337 338 return result; 339 } 340 341 /** 342 * Shallow stress balance class that reads `u` and `v` fields from a 343 * file and holds them constant. 344 * 345 * The only use I can think of right now is testing. 346 */ 347 PrescribedSliding::PrescribedSliding(IceGrid::ConstPtr g) 348 : ZeroSliding(g) { 349 // empty 350 } 351 352 PrescribedSliding::~PrescribedSliding() { 353 // empty 354 } 355 356 void PrescribedSliding::update(const Inputs &inputs, bool full_update) { 357 (void) inputs; 358 if (full_update) { 359 m_basal_frictional_heating.set(0.0); 360 } 361 } 362 363 void PrescribedSliding::init_impl() { 364 ShallowStressBalance::init_impl(); 365 366 auto input_filename = m_config->get_string("stress_balance.prescribed_sliding.file"); 367 368 if (input_filename.empty()) { 369 throw RuntimeError(PISM_ERROR_LOCATION, 370 "stress_balance.prescribed_sliding.file is required."); 371 } 372 373 m_velocity.regrid(input_filename, CRITICAL); 374 } 375 376 SSB_beta::SSB_beta(const ShallowStressBalance *m) 377 : Diag<ShallowStressBalance>(m) { 378 // set metadata: 379 m_vars = {SpatialVariableMetadata(m_sys, "beta")}; 380 381 set_attrs("basal drag coefficient", "", "Pa s / m", "Pa s / m", 0); 382 } 383 384 IceModelVec::Ptr SSB_beta::compute_impl() const { 385 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "beta", WITHOUT_GHOSTS)); 386 result->metadata(0) = m_vars[0]; 387 388 const IceModelVec2S *tauc = m_grid->variables().get_2d_scalar("tauc"); 389 390 const IceBasalResistancePlasticLaw *basal_sliding_law = model->sliding_law(); 391 392 const IceModelVec2V &velocity = model->velocity(); 393 394 IceModelVec::AccessList list{tauc, &velocity, result.get()}; 395 for (Points p(*m_grid); p; p.next()) { 396 const int i = p.i(), j = p.j(); 397 398 (*result)(i,j) = basal_sliding_law->drag((*tauc)(i,j), velocity(i,j).u, velocity(i,j).v); 399 } 400 401 return result; 402 } 403 404 } // end of namespace stressbalance 405 } // end of namespace pism