output_extra.cc (13349B)
1 /* Copyright (C) 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 <netcdf.h> 21 #ifdef NC_HAVE_META_H 22 #include <netcdf_meta.h> 23 #endif 24 25 #include "IceModel.hh" 26 27 #include "pism/util/pism_options.hh" 28 #include "pism/util/pism_utilities.hh" 29 #include "pism/util/Profiling.hh" 30 31 namespace pism { 32 33 //! Computes the maximum time-step we can take and still hit all `-extra_times`. 34 MaxTimestep IceModel::extras_max_timestep(double my_t) { 35 36 if ((not m_save_extra) or 37 (not m_config->get_flag("time_stepping.hit_extra_times"))) { 38 return MaxTimestep("reporting (-extra_times)"); 39 } 40 41 return reporting_max_timestep(m_extra_times, my_t, "reporting (-extra_times)"); 42 } 43 44 static std::set<std::string> process_extra_shortcuts(const Config &config, 45 const std::set<std::string> &input) { 46 std::set<std::string> result = input; 47 48 // process shortcuts 49 if (result.find("amount_fluxes") != result.end()) { 50 result.erase("amount_fluxes"); 51 result.insert("tendency_of_ice_amount"); 52 result.insert("tendency_of_ice_amount_due_to_basal_mass_flux"); 53 result.insert("tendency_of_ice_amount_due_to_conservation_error"); 54 result.insert("tendency_of_ice_amount_due_to_discharge"); 55 result.insert("tendency_of_ice_amount_due_to_flow"); 56 result.insert("tendency_of_ice_amount_due_to_surface_mass_flux"); 57 } 58 59 if (result.find("mass_fluxes") != result.end()) { 60 result.erase("mass_fluxes"); 61 result.insert("tendency_of_ice_mass"); 62 result.insert("tendency_of_ice_mass_due_to_basal_mass_flux"); 63 result.insert("tendency_of_ice_mass_due_to_conservation_error"); 64 result.insert("tendency_of_ice_mass_due_to_discharge"); 65 result.insert("tendency_of_ice_mass_due_to_flow"); 66 result.insert("tendency_of_ice_mass_due_to_surface_mass_flux"); 67 } 68 69 if (result.find("pdd_fluxes") != result.end()) { 70 result.erase("pdd_fluxes"); 71 result.insert("surface_accumulation_flux"); 72 result.insert("surface_runoff_flux"); 73 result.insert("surface_melt_flux"); 74 } 75 76 if (result.find("pdd_rates") != result.end()) { 77 result.erase("pdd_rates"); 78 result.insert("surface_accumulation_rate"); 79 result.insert("surface_runoff_rate"); 80 result.insert("surface_melt_rate"); 81 } 82 83 if (result.find("hydrology_fluxes") != result.end()) { 84 result.erase("hydrology_fluxes"); 85 result.insert("tendency_of_subglacial_water_mass"); 86 result.insert("tendency_of_subglacial_water_mass_due_to_input"); 87 result.insert("tendency_of_subglacial_water_mass_due_to_flow"); 88 result.insert("tendency_of_subglacial_water_mass_due_to_conservation_error"); 89 result.insert("tendency_of_subglacial_water_mass_at_grounded_margins"); 90 result.insert("tendency_of_subglacial_water_mass_at_grounding_line"); 91 result.insert("tendency_of_subglacial_water_mass_at_domain_boundary"); 92 } 93 94 if (result.find("ismip6") != result.end()) { 95 96 const char *flag_name = "output.ISMIP6"; 97 98 if (not config.get_flag(flag_name)) { 99 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Please set %s to save ISMIP6 diagnostics " 100 "(-extra_vars ismip6).", flag_name); 101 } 102 103 result.erase("ismip6"); 104 for (auto v : set_split(config.get_string("output.ISMIP6_extra_variables"), ',')) { 105 result.insert(v); 106 } 107 } 108 109 return result; 110 } 111 112 //! Initialize the code saving spatially-variable diagnostic quantities. 113 void IceModel::init_extras() { 114 115 m_last_extra = 0; // will be set in write_extras() 116 m_next_extra = 0; 117 118 m_extra_filename = m_config->get_string("output.extra.file"); 119 std::string times = m_config->get_string("output.extra.times"); 120 std::string vars = m_config->get_string("output.extra.vars"); 121 bool split = m_config->get_flag("output.extra.split"); 122 bool append = m_config->get_flag("output.extra.append"); 123 124 bool extra_file_set = not m_extra_filename.empty(); 125 bool times_set = not times.empty(); 126 127 if (extra_file_set ^ times_set) { 128 throw RuntimeError(PISM_ERROR_LOCATION, 129 "you need to set both output.extra.file and output.extra.times" 130 " to save spatial time-series."); 131 } 132 133 if (not extra_file_set and not times_set) { 134 m_save_extra = false; 135 return; 136 } 137 138 try { 139 m_extra_times = m_time->parse_times(times); 140 } catch (RuntimeError &e) { 141 e.add_context("parsing the output.extra.times argument %s", times.c_str()); 142 throw; 143 } 144 145 if (m_extra_times.size() == 0) { 146 throw RuntimeError(PISM_ERROR_LOCATION, "output.extra.times cannot be empty"); 147 } 148 149 if (append and split) { 150 throw RuntimeError(PISM_ERROR_LOCATION, 151 "both output.extra.split and output.extra.append are set."); 152 } 153 154 if (append) { 155 File file(m_grid->com, m_extra_filename, PISM_NETCDF3, PISM_READONLY); 156 157 std::string time_name = m_config->get_string("time.dimension_name"); 158 if (file.find_variable(time_name)) { 159 double time_max = vector_max(file.read_dimension(time_name)); 160 161 while (m_next_extra + 1 < m_extra_times.size() && m_extra_times[m_next_extra + 1] < time_max) { 162 m_next_extra++; 163 } 164 165 if (m_next_extra > 0) { 166 m_log->message(2, 167 "skipping times before the last record in %s (at %s)\n", 168 m_extra_filename.c_str(), m_time->date(time_max).c_str()); 169 } 170 171 // discard requested times before the beginning of the run 172 std::vector<double> tmp(m_extra_times.size() - m_next_extra); 173 for (unsigned int k = 0; k < tmp.size(); ++k) { 174 tmp[k] = m_extra_times[m_next_extra + k]; 175 } 176 177 m_extra_times = tmp; 178 m_next_extra = 0; 179 } 180 } 181 182 m_save_extra = true; 183 m_extra_file_is_ready = false; 184 m_split_extra = false; 185 186 if (split) { 187 m_split_extra = true; 188 m_log->message(2, "saving spatial time-series to '%s+year.nc'; ", 189 m_extra_filename.c_str()); 190 } else { 191 if (not ends_with(m_extra_filename, ".nc")) { 192 m_log->message(2, 193 "PISM WARNING: spatial time-series file name '%s' does not have the '.nc' suffix!\n", 194 m_extra_filename.c_str()); 195 } 196 m_log->message(2, "saving spatial time-series to '%s'; ", 197 m_extra_filename.c_str()); 198 } 199 200 m_log->message(2, "times requested: %s\n", times.c_str()); 201 202 if (m_extra_times.size() > 500) { 203 m_log->message(2, 204 "PISM WARNING: more than 500 times requested. This might fill your hard-drive!\n"); 205 } 206 207 #ifdef NC_HAVE_META_H 208 { 209 if (100 * NC_VERSION_MAJOR + 10 * NC_VERSION_MINOR + NC_VERSION_PATCH < 473) { 210 if (m_extra_times.size() > 5000 and m_config->get_string("output.format") == "netcdf4_parallel") { 211 throw RuntimeError(PISM_ERROR_LOCATION, 212 "more than 5000 times requested." 213 "Please use -extra_split to avoid a crash caused by a bug in NetCDF versions older than 4.7.3.\n" 214 "Alternatively\n" 215 "- split this simulation into several runs and then concatenate results\n" 216 "- select a different output.format value\n" 217 "- upgrade NetCDF to 4.7.3"); 218 } 219 } 220 } 221 #endif 222 223 if (not vars.empty()) { 224 m_extra_vars = process_extra_shortcuts(*m_config, set_split(vars, ',')); 225 m_log->message(2, "variables requested: %s\n", vars.c_str()); 226 } else { 227 m_log->message(2, 228 "PISM WARNING: output.extra.vars was not set. Writing the model state...\n"); 229 } // end of the else clause after "if (extra_vars_set)" 230 } 231 232 //! Write spatially-variable diagnostic quantities. 233 void IceModel::write_extras() { 234 double saving_after = -1.0e30; // initialize to avoid compiler warning; this 235 // value is never used, because saving_after 236 // is only used if save_now == true, and in 237 // this case saving_after is guaranteed to be 238 // initialized. See the code below. 239 char filename[PETSC_MAX_PATH_LEN]; 240 unsigned int current_extra; 241 // determine if the user set the -save_at and -save_to options 242 if (not m_save_extra) { 243 return; 244 } 245 246 double current_time = m_time->current(); 247 248 // do we need to save *now*? 249 if (m_next_extra < m_extra_times.size() and 250 (current_time >= m_extra_times[m_next_extra] or 251 fabs(current_time - m_extra_times[m_next_extra]) < 1.0)) { 252 // the condition above is "true" if we passed a requested time or got to 253 // within 1 second from it 254 255 current_extra = m_next_extra; 256 257 // update next_extra 258 while (m_next_extra < m_extra_times.size() and 259 (m_extra_times[m_next_extra] <= current_time or 260 fabs(current_time - m_extra_times[m_next_extra]) < 1.0)) { 261 m_next_extra++; 262 } 263 264 saving_after = m_extra_times[current_extra]; 265 } else { 266 return; 267 } 268 269 if (current_extra == 0) { 270 // The first time defines the left end-point of the first reporting interval; we don't write a 271 // report at this time. 272 273 // Re-initialize last_extra (the correct value is not known at the time init_extras() is 274 // called). 275 m_last_extra = current_time; 276 277 // ISMIP6 runs need to save diagnostics at the beginning of the run 278 if (not m_config->get_flag("output.ISMIP6")) { 279 return; 280 } 281 } 282 283 if (saving_after < m_time->start()) { 284 // Suppose a user tells PISM to write data at times 0:1000:10000. Suppose 285 // also that PISM writes a backup file at year 2500 and gets stopped. 286 // 287 // When restarted, PISM will decide that it's time to write data for time 288 // 2000, but 289 // * that record was written already and 290 // * PISM will end up writing at year 2500, producing a file containing one 291 // more record than necessary. 292 // 293 // This check makes sure that this never happens. 294 return; 295 } 296 297 if (m_split_extra) { 298 m_extra_file_is_ready = false; // each time-series record is written to a separate file 299 snprintf(filename, PETSC_MAX_PATH_LEN, "%s_%s.nc", 300 m_extra_filename.c_str(), m_time->date().c_str()); 301 } else { 302 strncpy(filename, m_extra_filename.c_str(), PETSC_MAX_PATH_LEN); 303 } 304 305 m_log->message(3, 306 "saving spatial time-series to %s at %s\n", 307 filename, m_time->date().c_str()); 308 309 // default behavior is to move the file aside if it exists already; option allows appending 310 bool append = m_config->get_flag("output.extra.append"); 311 IO_Mode mode = m_extra_file_is_ready or append ? PISM_READWRITE : PISM_READWRITE_MOVE; 312 313 const Profiling &profiling = m_ctx->profiling(); 314 profiling.begin("io.extra_file"); 315 { 316 if (not m_extra_file) { 317 m_extra_file.reset(new File(m_grid->com, 318 filename, 319 string_to_backend(m_config->get_string("output.format")), 320 mode, 321 m_ctx->pio_iosys_id())); 322 } 323 324 std::string time_name = m_config->get_string("time.dimension_name"); 325 326 if (not m_extra_file_is_ready) { 327 // Prepare the file: 328 io::define_time(*m_extra_file, *m_ctx); 329 m_extra_file->write_attribute(time_name, "bounds", "time_bounds"); 330 331 io::define_time_bounds(m_extra_bounds, *m_extra_file); 332 333 write_metadata(*m_extra_file, WRITE_MAPPING, PREPEND_HISTORY); 334 335 m_extra_file_is_ready = true; 336 } 337 338 write_run_stats(*m_extra_file); 339 340 save_variables(*m_extra_file, 341 m_extra_vars.empty() ? INCLUDE_MODEL_STATE : JUST_DIAGNOSTICS, 342 m_extra_vars, 343 0.5 * (m_last_extra + current_time), // use the mid-point of the 344 // current reporting interval 345 PISM_FLOAT); 346 347 // Get the length of the time dimension *after* it is appended to. 348 unsigned int time_length = m_extra_file->dimension_length(time_name); 349 size_t time_start = time_length > 0 ? static_cast<size_t>(time_length - 1) : 0; 350 351 io::write_time_bounds(*m_extra_file, m_extra_bounds, 352 time_start, {m_last_extra, current_time}); 353 // make sure all changes are written 354 m_extra_file->sync(); 355 } 356 profiling.end("io.extra_file"); 357 358 flush_timeseries(); 359 360 if (m_split_extra) { 361 // each record is saved to a new file, so we can close this one 362 m_extra_file.reset(nullptr); 363 } 364 365 m_last_extra = current_time; 366 367 // reset accumulators in diagnostics that compute time averaged quantities 368 reset_diagnostics(); 369 } 370 371 } // end of namespace pism