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