toutput_extra.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
toutput_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