EnergyModel.cc (13772B)
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 20 #include "EnergyModel.hh" 21 #include "pism/util/MaxTimestep.hh" 22 #include "pism/stressbalance/StressBalance.hh" 23 #include "pism/util/io/File.hh" 24 #include "pism/util/Vars.hh" 25 #include "utilities.hh" 26 #include "pism/util/EnthalpyConverter.hh" 27 #include "bootstrapping.hh" 28 #include "pism/util/pism_utilities.hh" 29 #include "pism/util/error_handling.hh" 30 #include "pism/util/IceModelVec2CellType.hh" 31 #include "pism/util/pism_options.hh" 32 #include "pism/util/Profiling.hh" 33 34 namespace pism { 35 namespace energy { 36 37 static void check_input(const IceModelVec *ptr, const char *name) { 38 if (ptr == NULL) { 39 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "energy balance model input %s was not provided", name); 40 } 41 } 42 43 Inputs::Inputs() { 44 basal_frictional_heating = NULL; 45 basal_heat_flux = NULL; 46 cell_type = NULL; 47 ice_thickness = NULL; 48 surface_liquid_fraction = NULL; 49 shelf_base_temp = NULL; 50 surface_temp = NULL; 51 till_water_thickness = NULL; 52 53 volumetric_heating_rate = NULL; 54 u3 = NULL; 55 v3 = NULL; 56 w3 = NULL; 57 58 no_model_mask = NULL; 59 } 60 61 void Inputs::check() const { 62 check_input(cell_type, "cell_type"); 63 check_input(basal_frictional_heating, "basal_frictional_heating"); 64 check_input(basal_heat_flux, "basal_heat_flux"); 65 check_input(ice_thickness, "ice_thickness"); 66 check_input(surface_liquid_fraction, "surface_liquid_fraction"); 67 check_input(shelf_base_temp, "shelf_base_temp"); 68 check_input(surface_temp, "surface_temp"); 69 check_input(till_water_thickness, "till_water_thickness"); 70 71 check_input(volumetric_heating_rate, "volumetric_heating_rate"); 72 check_input(u3, "u3"); 73 check_input(v3, "v3"); 74 check_input(w3, "w3"); 75 } 76 77 EnergyModelStats::EnergyModelStats() { 78 bulge_counter = 0; 79 reduced_accuracy_counter = 0; 80 low_temperature_counter = 0; 81 liquified_ice_volume = 0.0; 82 } 83 84 EnergyModelStats& EnergyModelStats::operator+=(const EnergyModelStats &other) { 85 bulge_counter += other.bulge_counter; 86 reduced_accuracy_counter += other.reduced_accuracy_counter; 87 low_temperature_counter += other.low_temperature_counter; 88 liquified_ice_volume += other.liquified_ice_volume; 89 return *this; 90 } 91 92 93 bool marginal(const IceModelVec2S &thickness, int i, int j, double threshold) { 94 int 95 n = j + 1, 96 e = i + 1, 97 s = j - 1, 98 w = i - 1; 99 100 const double 101 N = thickness(i, n), 102 E = thickness(e, j), 103 S = thickness(i, s), 104 W = thickness(w, j), 105 NW = thickness(w, n), 106 SW = thickness(w, s), 107 NE = thickness(e, n), 108 SE = thickness(e, s); 109 110 return ((E < threshold) or 111 (NE < threshold) or 112 (N < threshold) or 113 (NW < threshold) or 114 (W < threshold) or 115 (SW < threshold) or 116 (S < threshold) or 117 (SE < threshold)); 118 } 119 120 121 void EnergyModelStats::sum(MPI_Comm com) { 122 bulge_counter = GlobalSum(com, bulge_counter); 123 reduced_accuracy_counter = GlobalSum(com, reduced_accuracy_counter); 124 low_temperature_counter = GlobalSum(com, low_temperature_counter); 125 liquified_ice_volume = GlobalSum(com, liquified_ice_volume); 126 } 127 128 129 EnergyModel::EnergyModel(IceGrid::ConstPtr grid, 130 stressbalance::StressBalance *stress_balance) 131 : Component(grid), m_stress_balance(stress_balance) { 132 133 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width"); 134 135 { 136 m_ice_enthalpy.create(m_grid, "enthalpy", WITH_GHOSTS, WIDE_STENCIL); 137 // POSSIBLE standard name = land_ice_enthalpy 138 m_ice_enthalpy.set_attrs("model_state", 139 "ice enthalpy (includes sensible heat, latent heat, pressure)", 140 "J kg-1", "J kg-1", "", 0); 141 } 142 143 { 144 m_basal_melt_rate.create(m_grid, "basal_melt_rate_grounded", WITHOUT_GHOSTS); 145 // ghosted to allow the "redundant" computation of tauc 146 m_basal_melt_rate.set_attrs("model_state", 147 "ice basal melt rate from energy conservation, in ice thickness per time (valid in grounded areas)", 148 "m s-1", "m year-1", "", 0); 149 // We could use land_ice_basal_melt_rate, but that way both basal_melt_rate_grounded and bmelt 150 // have this standard name. 151 m_basal_melt_rate.metadata().set_string("comment", "positive basal melt rate corresponds to ice loss"); 152 } 153 154 // a 3d work vector 155 { 156 m_work.create(m_grid, "work_vector", WITHOUT_GHOSTS); 157 m_work.set_attrs("internal", 158 "usually new values of temperature or enthalpy during time step", 159 "", "", "", 0); 160 } 161 } 162 163 void EnergyModel::init_enthalpy(const File &input_file, bool do_regrid, int record) { 164 165 if (input_file.find_variable("enthalpy")) { 166 if (do_regrid) { 167 m_ice_enthalpy.regrid(input_file, CRITICAL); 168 } else { 169 m_ice_enthalpy.read(input_file, record); 170 } 171 } else if (input_file.find_variable("temp")) { 172 IceModelVec3 173 &temp = m_work, 174 &liqfrac = m_ice_enthalpy; 175 176 { 177 temp.set_name("temp"); 178 temp.metadata(0).set_name("temp"); 179 temp.set_attrs("temporary", "ice temperature", 180 "Kelvin", "Kelvin", "land_ice_temperature", 0); 181 182 if (do_regrid) { 183 temp.regrid(input_file, CRITICAL); 184 } else { 185 temp.read(input_file, record); 186 } 187 } 188 189 const IceModelVec2S & ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness"); 190 191 if (input_file.find_variable("liqfrac")) { 192 SpatialVariableMetadata enthalpy_metadata = m_ice_enthalpy.metadata(); 193 194 liqfrac.set_name("liqfrac"); 195 liqfrac.metadata(0).set_name("liqfrac"); 196 liqfrac.set_attrs("temporary", "ice liquid water fraction", 197 "1", "1", "", 0); 198 199 if (do_regrid) { 200 liqfrac.regrid(input_file, CRITICAL); 201 } else { 202 liqfrac.read(input_file, record); 203 } 204 205 m_ice_enthalpy.set_name(enthalpy_metadata.get_name()); 206 m_ice_enthalpy.metadata() = enthalpy_metadata; 207 208 m_log->message(2, 209 " - Computing enthalpy using ice temperature," 210 " liquid water fraction and thickness...\n"); 211 compute_enthalpy(temp, liqfrac, ice_thickness, m_ice_enthalpy); 212 } else { 213 m_log->message(2, 214 " - Computing enthalpy using ice temperature and thickness " 215 "and assuming zero liquid water fraction...\n"); 216 compute_enthalpy_cold(temp, ice_thickness, m_ice_enthalpy); 217 } 218 } else { 219 throw RuntimeError::formatted(PISM_ERROR_LOCATION, 220 "neither enthalpy nor temperature was found in '%s'.\n", 221 input_file.filename().c_str()); 222 } 223 } 224 225 /*! 226 * The `-regrid_file` may contain enthalpy, temperature, or *both* temperature and liquid water 227 * fraction. 228 */ 229 void EnergyModel::regrid_enthalpy() { 230 231 auto regrid_filename = m_config->get_string("input.regrid.file"); 232 auto regrid_vars = set_split(m_config->get_string("input.regrid.vars"), ','); 233 234 235 if (regrid_filename.empty()) { 236 return; 237 } 238 239 std::string enthalpy_name = m_ice_enthalpy.metadata().get_name(); 240 241 if (regrid_vars.empty() or member(enthalpy_name, regrid_vars)) { 242 File regrid_file(m_grid->com, regrid_filename, PISM_GUESS, PISM_READONLY); 243 init_enthalpy(regrid_file, true, 0); 244 } 245 } 246 247 248 void EnergyModel::restart(const File &input_file, int record) { 249 this->restart_impl(input_file, record); 250 } 251 252 void EnergyModel::bootstrap(const File &input_file, 253 const IceModelVec2S &ice_thickness, 254 const IceModelVec2S &surface_temperature, 255 const IceModelVec2S &climatic_mass_balance, 256 const IceModelVec2S &basal_heat_flux) { 257 this->bootstrap_impl(input_file, 258 ice_thickness, surface_temperature, 259 climatic_mass_balance, basal_heat_flux); 260 } 261 262 void EnergyModel::initialize(const IceModelVec2S &basal_melt_rate, 263 const IceModelVec2S &ice_thickness, 264 const IceModelVec2S &surface_temperature, 265 const IceModelVec2S &climatic_mass_balance, 266 const IceModelVec2S &basal_heat_flux) { 267 this->initialize_impl(basal_melt_rate, 268 ice_thickness, 269 surface_temperature, 270 climatic_mass_balance, 271 basal_heat_flux); 272 } 273 274 void EnergyModel::update(double t, double dt, const Inputs &inputs) { 275 // reset standard out flags at the beginning of every time step 276 m_stdout_flags = ""; 277 m_stats = EnergyModelStats(); 278 279 const Profiling &profiling = m_grid->ctx()->profiling(); 280 281 profiling.begin("ice_energy"); 282 { 283 // this call should fill m_work with new values of enthalpy 284 this->update_impl(t, dt, inputs); 285 286 m_work.update_ghosts(m_ice_enthalpy); 287 } 288 profiling.end("ice_energy"); 289 290 // globalize m_stats and update m_stdout_flags 291 { 292 char buffer[50] = ""; 293 294 m_stats.sum(m_grid->com); 295 296 if (m_stats.reduced_accuracy_counter > 0.0) { // count of when BOMBPROOF switches to lower accuracy 297 const double reduced_accuracy_percentage = 100.0 * (m_stats.reduced_accuracy_counter / (m_grid->Mx() * m_grid->My())); 298 const double reporting_threshold = 5.0; // only report if above 5% 299 300 if (reduced_accuracy_percentage > reporting_threshold and m_log->get_threshold() > 2) { 301 snprintf(buffer, 50, " [BPsacr=%.4f%%] ", reduced_accuracy_percentage); 302 m_stdout_flags = buffer + m_stdout_flags; 303 } 304 } 305 306 if (m_stats.bulge_counter > 0) { 307 // count of when advection bulges are limited; frequently it is identically zero 308 snprintf(buffer, 50, " BULGE=%d ", m_stats.bulge_counter); 309 m_stdout_flags = buffer + m_stdout_flags; 310 } 311 } 312 } 313 314 MaxTimestep EnergyModel::max_timestep_impl(double t) const { 315 // silence a compiler warning 316 (void) t; 317 318 if (m_stress_balance == NULL) { 319 throw RuntimeError::formatted(PISM_ERROR_LOCATION, 320 "EnergyModel: no stress balance provided." 321 " Cannot compute max. time step."); 322 } 323 324 return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "energy"); 325 } 326 327 const std::string& EnergyModel::stdout_flags() const { 328 return m_stdout_flags; 329 } 330 331 const EnergyModelStats& EnergyModel::stats() const { 332 return m_stats; 333 } 334 335 const IceModelVec3 & EnergyModel::enthalpy() const { 336 return m_ice_enthalpy; 337 } 338 339 /*! @brief Basal melt rate in grounded areas. (It is set to zero elsewhere.) */ 340 const IceModelVec2S & EnergyModel::basal_melt_rate() const { 341 return m_basal_melt_rate; 342 } 343 344 /*! @brief Ice loss "flux" due to ice liquefaction. */ 345 class LiquifiedIceFlux : public TSDiag<TSFluxDiagnostic,EnergyModel> { 346 public: 347 LiquifiedIceFlux(const EnergyModel *m) 348 : TSDiag<TSFluxDiagnostic, EnergyModel>(m, "liquified_ice_flux") { 349 350 set_units("m3 / second", "m3 / year"); 351 m_ts.variable().set_string("long_name", 352 "rate of ice loss due to liquefaction," 353 " averaged over the reporting interval"); 354 m_ts.variable().set_string("comment", "positive means ice loss"); 355 m_ts.variable().set_string("cell_methods", "time: mean"); 356 } 357 protected: 358 double compute() { 359 // liquified ice volume during the last time step 360 return model->stats().liquified_ice_volume; 361 } 362 }; 363 364 namespace diagnostics { 365 /*! @brief Report ice enthalpy. */ 366 class Enthalpy : public Diag<EnergyModel> 367 { 368 public: 369 Enthalpy(const EnergyModel *m) 370 : Diag<EnergyModel>(m) { 371 m_vars = {model->enthalpy().metadata()}; 372 } 373 374 protected: 375 IceModelVec::Ptr compute_impl() const { 376 377 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "enthalpy", WITHOUT_GHOSTS)); 378 result->metadata(0) = m_vars[0]; 379 380 const IceModelVec3 &input = model->enthalpy(); 381 382 // FIXME: implement IceModelVec3::copy_from() 383 384 IceModelVec::AccessList list {result.get(), &input}; 385 ParallelSection loop(m_grid->com); 386 try { 387 for (Points p(*m_grid); p; p.next()) { 388 const int i = p.i(), j = p.j(); 389 390 result->set_column(i, j, input.get_column(i, j)); 391 } 392 } catch (...) { 393 loop.failed(); 394 } 395 loop.check(); 396 397 398 return result; 399 } 400 }; 401 402 } // end of namespace diagnostics 403 404 DiagnosticList EnergyModel::diagnostics_impl() const { 405 DiagnosticList result; 406 result = { 407 {"enthalpy", Diagnostic::Ptr(new diagnostics::Enthalpy(this))}, 408 {"basal_melt_rate_grounded", Diagnostic::wrap(m_basal_melt_rate)} 409 }; 410 return result; 411 } 412 413 TSDiagnosticList EnergyModel::ts_diagnostics_impl() const { 414 return { 415 {"liquified_ice_flux", TSDiagnostic::Ptr(new LiquifiedIceFlux(this))} 416 }; 417 } 418 419 } // end of namespace energy 420 421 } // end of namespace pism