PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
output_extra.cc
Go to the documentation of this file.
1 /* Copyright (C) 2017, 2018, 2019, 2020, 2021, 2023 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 "pism/icemodel/IceModel.hh"
26 
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/util/Profiling.hh"
29 
30 namespace pism {
31 
32 //! Computes the maximum time-step we can take and still hit all `-extra_times`.
34 
35  if ((not m_save_extra) or
36  (not m_config->get_flag("time_stepping.hit_extra_times"))) {
37  return MaxTimestep("reporting (-extra_times)");
38  }
39 
40  double eps = m_config->get_number("time_stepping.resolution");
41 
42  return reporting_max_timestep(m_extra_times, my_t, eps,
43  "reporting (-extra_times)");
44 }
45 
46 static std::set<std::string> process_extra_shortcuts(const Config &config,
47  const std::set<std::string> &input) {
48  std::set<std::string> result = input;
49 
50  // process shortcuts
51  if (result.find("amount_fluxes") != result.end()) {
52  result.erase("amount_fluxes");
53  result.insert("tendency_of_ice_amount");
54  result.insert("tendency_of_ice_amount_due_to_basal_mass_flux");
55  result.insert("tendency_of_ice_amount_due_to_conservation_error");
56  result.insert("tendency_of_ice_amount_due_to_discharge");
57  result.insert("tendency_of_ice_amount_due_to_flow");
58  result.insert("tendency_of_ice_amount_due_to_surface_mass_flux");
59  }
60 
61  if (result.find("mass_fluxes") != result.end()) {
62  result.erase("mass_fluxes");
63  result.insert("tendency_of_ice_mass");
64  result.insert("tendency_of_ice_mass_due_to_basal_mass_flux");
65  result.insert("tendency_of_ice_mass_due_to_conservation_error");
66  result.insert("tendency_of_ice_mass_due_to_discharge");
67  result.insert("tendency_of_ice_mass_due_to_flow");
68  result.insert("tendency_of_ice_mass_due_to_surface_mass_flux");
69  }
70 
71  if (result.find("pdd_fluxes") != result.end()) {
72  result.erase("pdd_fluxes");
73  result.insert("surface_accumulation_flux");
74  result.insert("surface_runoff_flux");
75  result.insert("surface_melt_flux");
76  }
77 
78  if (result.find("pdd_rates") != result.end()) {
79  result.erase("pdd_rates");
80  result.insert("surface_accumulation_rate");
81  result.insert("surface_runoff_rate");
82  result.insert("surface_melt_rate");
83  }
84 
85  if (result.find("hydrology_fluxes") != result.end()) {
86  result.erase("hydrology_fluxes");
87  result.insert("tendency_of_subglacial_water_mass");
88  result.insert("tendency_of_subglacial_water_mass_due_to_input");
89  result.insert("tendency_of_subglacial_water_mass_due_to_flow");
90  result.insert("tendency_of_subglacial_water_mass_due_to_conservation_error");
91  result.insert("tendency_of_subglacial_water_mass_at_grounded_margins");
92  result.insert("tendency_of_subglacial_water_mass_at_grounding_line");
93  result.insert("tendency_of_subglacial_water_mass_at_domain_boundary");
94  }
95 
96  if (result.find("ismip6") != result.end()) {
97 
98  const char *flag_name = "output.ISMIP6";
99 
100  if (not config.get_flag(flag_name)) {
101  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Please set %s to save ISMIP6 diagnostics "
102  "(-extra_vars ismip6).", flag_name);
103  }
104 
105  result.erase("ismip6");
106  for (const auto& v : set_split(config.get_string("output.ISMIP6_extra_variables"), ',')) {
107  result.insert(v);
108  }
109  }
110 
111  return result;
112 }
113 
114 //! Initialize the code saving spatially-variable diagnostic quantities.
116 
117  m_last_extra = 0; // will be set in write_extras()
118  m_next_extra = 0;
119 
120  m_extra_filename = m_config->get_string("output.extra.file");
121  std::string times = m_config->get_string("output.extra.times");
122  std::string vars = m_config->get_string("output.extra.vars");
123  bool split = m_config->get_flag("output.extra.split");
124  bool append = m_config->get_flag("output.extra.append");
125 
126  bool extra_file_set = not m_extra_filename.empty();
127  bool times_set = not times.empty();
128 
129  if (extra_file_set ^ times_set) {
131  "you need to set both output.extra.file and output.extra.times"
132  " to save spatial time-series.");
133  }
134 
135  if (not extra_file_set and not times_set) {
136  m_save_extra = false;
137  return;
138  }
139 
140  try {
141  m_extra_times = m_time->parse_times(times);
142  } catch (RuntimeError &e) {
143  e.add_context("parsing the output.extra.times argument %s", times.c_str());
144  throw;
145  }
146 
147  if (m_extra_times.size() == 0) {
148  throw RuntimeError(PISM_ERROR_LOCATION, "output.extra.times cannot be empty");
149  }
150 
151  if (append and split) {
153  "both output.extra.split and output.extra.append are set.");
154  }
155 
156  if (append) {
158 
159  std::string time_name = m_config->get_string("time.dimension_name");
160  if (file.find_variable(time_name)) {
161  double time_max = vector_max(file.read_dimension(time_name));
162 
163  while (m_next_extra + 1 < m_extra_times.size() && m_extra_times[m_next_extra + 1] < time_max) {
164  m_next_extra++;
165  }
166 
167  if (m_next_extra > 0) {
168  m_log->message(2,
169  "skipping times before the last record in %s (at %s)\n",
170  m_extra_filename.c_str(), m_time->date(time_max).c_str());
171  }
172 
173  // discard requested times before the beginning of the run
174  std::vector<double> tmp(m_extra_times.size() - m_next_extra);
175  for (unsigned int k = 0; k < tmp.size(); ++k) {
176  tmp[k] = m_extra_times[m_next_extra + k];
177  }
178 
179  m_extra_times = tmp;
180  m_next_extra = 0;
181  }
182  }
183 
184  m_save_extra = true;
185  m_extra_file_is_ready = false;
186  m_split_extra = false;
187 
188  if (split) {
189  m_split_extra = true;
190  m_log->message(2, "saving spatial time-series to '%s+year.nc'; ",
191  m_extra_filename.c_str());
192  } else {
193  if (not ends_with(m_extra_filename, ".nc")) {
194  m_log->message(2,
195  "PISM WARNING: spatial time-series file name '%s' does not have the '.nc' suffix!\n",
196  m_extra_filename.c_str());
197  }
198  m_log->message(2, "saving spatial time-series to '%s'; ",
199  m_extra_filename.c_str());
200  }
201 
202  m_log->message(2, "times requested: %s\n", times.c_str());
203 
204  if (m_extra_times.size() > 500) {
205  m_log->message(2,
206  "PISM WARNING: more than 500 times requested. This might fill your hard-drive!\n");
207  }
208 
209 #ifdef NC_HAVE_META_H
210  {
211  if (100 * NC_VERSION_MAJOR + 10 * NC_VERSION_MINOR + NC_VERSION_PATCH < 473) {
212  if (m_extra_times.size() > 5000 and m_config->get_string("output.format") == "netcdf4_parallel") {
214  "more than 5000 times requested."
215  "Please use -extra_split to avoid a crash caused by a bug in NetCDF versions older than 4.7.3.\n"
216  "Alternatively\n"
217  "- split this simulation into several runs and then concatenate results\n"
218  "- select a different output.format value\n"
219  "- upgrade NetCDF to 4.7.3");
220  }
221  }
222  }
223 #endif
224 
225  if (not vars.empty()) {
227  m_log->message(2, "variables requested: %s\n", vars.c_str());
228  } else {
229  m_log->message(2,
230  "PISM WARNING: output.extra.vars was not set. Writing the model state...\n");
231  } // end of the else clause after "if (extra_vars_set)"
232 }
233 
234 //! Write spatially-variable diagnostic quantities.
236  double saving_after = -1.0e30; // initialize to avoid compiler warning; this
237  // value is never used, because saving_after
238  // is only used if save_now == true, and in
239  // this case saving_after is guaranteed to be
240  // initialized. See the code below.
241  std::string filename;
242  unsigned int current_extra;
243  // determine if the user set the -save_at and -save_to options
244  if (not m_save_extra) {
245  return;
246  }
247 
248  const double time_resolution = m_config->get_number("time_stepping.resolution");
249  double current_time = m_time->current();
250 
251  // do we need to save *now*?
252  if (m_next_extra < m_extra_times.size() and
253  (current_time >= m_extra_times[m_next_extra] or
254  fabs(current_time - m_extra_times[m_next_extra]) < time_resolution)) {
255  // the condition above is "true" if we passed a requested time or got to
256  // within time_resolution seconds from it
257 
258  current_extra = m_next_extra;
259 
260  // update next_extra
261  while (m_next_extra < m_extra_times.size() and
262  (m_extra_times[m_next_extra] <= current_time or
263  fabs(current_time - m_extra_times[m_next_extra]) < time_resolution)) {
264  m_next_extra++;
265  }
266 
267  saving_after = m_extra_times[current_extra];
268  } else {
269  return;
270  }
271 
272  if (current_extra == 0) {
273  // The first time defines the left end-point of the first reporting interval; we don't write a
274  // report at this time.
275 
276  // Re-initialize last_extra (the correct value is not known at the time init_extras() is
277  // called).
278  m_last_extra = current_time;
279 
280  // ISMIP6 runs need to save diagnostics at the beginning of the run
281  if (not m_config->get_flag("output.ISMIP6")) {
282  return;
283  }
284  }
285 
286  if (saving_after < m_time->start()) {
287  // Suppose a user tells PISM to write data at times 0:1000:10000. Suppose
288  // also that PISM writes a backup file at year 2500 and gets stopped.
289  //
290  // When restarted, PISM will decide that it's time to write data for time
291  // 2000, but
292  // * that record was written already and
293  // * PISM will end up writing at year 2500, producing a file containing one
294  // more record than necessary.
295  //
296  // This check makes sure that this never happens.
297  return;
298  }
299 
300  if (m_split_extra) {
301  m_extra_file_is_ready = false; // each time-series record is written to a separate file
302  auto date_without_spaces = replace_character(m_time->date(m_time->current()), ' ', '_');
303  filename = pism::printf("%s_%s.nc",
304  m_extra_filename.c_str(),
305  date_without_spaces.c_str());
306  } else {
307  filename = m_extra_filename;
308  }
309 
310  m_log->message(3,
311  "saving spatial time-series to %s at %s\n",
312  filename.c_str(), m_time->date(m_time->current()).c_str());
313 
314  // default behavior is to move the file aside if it exists already; option allows appending
315  bool append = m_config->get_flag("output.extra.append");
317 
318  const Profiling &profiling = m_ctx->profiling();
319  profiling.begin("io.extra_file");
320  {
321  if (not m_extra_file) {
322  m_extra_file.reset(new File(m_grid->com,
323  filename,
324  string_to_backend(m_config->get_string("output.format")),
325  mode,
326  m_ctx->pio_iosys_id()));
327  }
328 
329  std::string time_name = m_config->get_string("time.dimension_name");
330 
331  if (not m_extra_file_is_ready) {
332  // Prepare the file:
334  m_extra_file->write_attribute(time_name, "bounds", "time_bounds");
335 
337  time_name, "nv", *m_extra_file, io::PISM_DOUBLE);
338 
340 
341  m_extra_file_is_ready = true;
342  }
343 
345 
348  m_extra_vars,
349  0.5 * (m_last_extra + current_time), // use the mid-point of the
350  // current reporting interval
352 
353  // Get the length of the time dimension *after* it is appended to.
354  unsigned int time_length = m_extra_file->dimension_length(time_name);
355  size_t time_start = time_length > 0 ? static_cast<size_t>(time_length - 1) : 0;
356 
358  time_start, {m_last_extra, current_time});
359  // make sure all changes are written
360  m_extra_file->sync();
361  }
362  profiling.end("io.extra_file");
363 
365 
366  if (m_split_extra) {
367  // each record is saved to a new file, so we can close this one
368  m_extra_file.reset(nullptr);
369  }
370 
371  m_last_extra = current_time;
372 
373  // reset accumulators in diagnostics that compute time averaged quantities
375 }
376 
377 } // end of namespace pism
std::string get_string(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
bool get_flag(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
Definition: File.cc:361
std::vector< double > read_dimension(const std::string &name) const
Get dimension data (a coordinate variable).
Definition: File.cc:603
High-level PISM I/O class.
Definition: File.hh:56
MaxTimestep extras_max_timestep(double my_t)
Computes the maximum time-step we can take and still hit all -extra_times.
Definition: output_extra.cc:33
VariableMetadata run_stats() const
Definition: utilities.cc:92
unsigned int m_next_extra
Definition: IceModel.hh:449
std::unique_ptr< File > m_extra_file
Definition: IceModel.hh:453
const Config::Ptr m_config
Configuration flags and parameters.
Definition: IceModel.hh:243
void write_extras()
Write spatially-variable diagnostic quantities.
bool m_save_extra
Definition: IceModel.hh:446
const Time::Ptr m_time
Time manager.
Definition: IceModel.hh:251
virtual void reset_diagnostics()
Definition: IceModel.cc:1063
void init_extras()
Initialize the code saving spatially-variable diagnostic quantities.
const std::shared_ptr< Context > m_ctx
Execution context.
Definition: IceModel.hh:245
const Logger::Ptr m_log
Logger.
Definition: IceModel.hh:249
void flush_timeseries()
Flush scalar time-series.
Definition: output_ts.cc:122
double m_last_extra
Definition: IceModel.hh:450
virtual void save_variables(const File &file, OutputKind kind, const std::set< std::string > &variables, double time, io::Type default_diagnostics_type=io::PISM_FLOAT) const
Definition: output.cc:169
std::set< std::string > m_extra_vars
Definition: IceModel.hh:451
virtual void write_metadata(const File &file, MappingTreatment mapping_flag, HistoryTreatment history_flag) const
Write time-independent metadata to a file.
Definition: output.cc:72
std::string m_extra_filename
Definition: IceModel.hh:447
bool m_extra_file_is_ready
Definition: IceModel.hh:446
std::vector< double > m_extra_times
Definition: IceModel.hh:448
bool m_split_extra
Definition: IceModel.hh:446
const std::shared_ptr< Grid > m_grid
Computational grid.
Definition: IceModel.hh:241
VariableMetadata m_extra_bounds
Definition: IceModel.hh:452
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
void begin(const char *name) const
Definition: Profiling.cc:75
void end(const char *name) const
Definition: Profiling.cc:91
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
#define PISM_ERROR_LOCATION
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READWRITE_MOVE
create a file for writing, move foo.nc to foo.nc~ if present
Definition: IO_Flags.hh:78
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
@ PISM_READWRITE
open an existing file for reading and writing
Definition: IO_Flags.hh:74
@ PISM_FLOAT
Definition: IO_Flags.hh:51
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
void define_time_bounds(const VariableMetadata &var, const std::string &dimension_name, const std::string &bounds_name, const File &file, io::Type nctype)
Definition: io_helpers.cc:951
void write_time_bounds(const File &file, const VariableMetadata &metadata, size_t t_start, const std::vector< double > &data)
Definition: io_helpers.cc:1057
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
Definition: io_helpers.cc:213
MaxTimestep reporting_max_timestep(const std::vector< double > &times, double t, double eps, const std::string &description)
Definition: output.cc:41
io::Backend string_to_backend(const std::string &backend)
Definition: File.cc:57
bool ends_with(const std::string &str, const std::string &suffix)
Returns true if str ends with suffix and false otherwise.
std::string printf(const char *format,...)
static const double k
Definition: exactTestP.cc:42
std::set< std::string > set_split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a set of strings.
static std::set< std::string > process_extra_shortcuts(const Config &config, const std::set< std::string > &input)
Definition: output_extra.cc:46
std::string replace_character(const std::string &input, char from, char to)
void write_run_stats(const File &file, const pism::VariableMetadata &stats)
Definition: output.cc:162
double vector_max(const std::vector< double > &input)
std::vector< std::string > split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a vector of strings.