PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Diagnostic.cc
Go to the documentation of this file.
1 /* Copyright (C) 2015, 2016, 2017, 2019, 2020, 2021, 2022 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 "Diagnostic.hh"
21 #include "pism/util/Time.hh"
22 #include "error_handling.hh"
23 #include "io/io_helpers.hh"
24 #include "pism/util/Logger.hh"
25 #include "pism/util/pism_utilities.hh"
26 #include "pism/util/Context.hh"
27 
28 namespace pism {
29 
31  : m_grid(grid),
32  m_sys(grid->ctx()->unit_system()),
33  m_config(grid->ctx()->config()),
34  m_fill_value(m_config->get_number("output.fill_value")) {
35  // empty
36 }
37 
38 void Diagnostic::update(double dt) {
39  this->update_impl(dt);
40 }
41 
42 void Diagnostic::update_impl(double dt) {
43  (void) dt;
44  // empty
45 }
46 
48  this->reset_impl();
49 }
50 
52  // empty
53 }
54 
55 /*!
56  * Convert from external (output) units to internal units.
57  */
58 double Diagnostic::to_internal(double x) const {
59  std::string
60  out = m_vars.at(0).get_string("glaciological_units");
61  std::string
62  in = m_vars.at(0).get_string("units");
63  return convert(m_sys, x, out, in);
64 }
65 
66 /*!
67  * Convert from internal to external (output) units.
68  */
69 double Diagnostic::to_external(double x) const {
70  std::string
71  out = m_vars.at(0).get_string("glaciological_units"),
72  in = m_vars.at(0).get_string("units");
73  return convert(m_sys, x, in, out);
74 }
75 
76 //! Get the number of NetCDF variables corresponding to a diagnostic quantity.
77 unsigned int Diagnostic::n_variables() const {
78  return m_vars.size();
79 }
80 
81 void Diagnostic::init(const File &input, unsigned int time) {
82  this->init_impl(input, time);
83 }
84 
85 void Diagnostic::define_state(const File &output) const {
86  this->define_state_impl(output);
87 }
88 
89 void Diagnostic::write_state(const File &output) const {
90  this->write_state_impl(output);
91 }
92 
93 void Diagnostic::init_impl(const File &input, unsigned int time) {
94  (void) input;
95  (void) time;
96  // empty
97 }
98 
99 void Diagnostic::define_state_impl(const File &output) const {
100  (void) output;
101  // empty
102 }
103 
104 void Diagnostic::write_state_impl(const File &output) const {
105  (void) output;
106  // empty
107 }
108 
109 //! Get a metadata object corresponding to variable number N.
111  if (N >= m_vars.size()) {
113  "variable metadata index %d is out of bounds",
114  N);
115  }
116 
117  return m_vars[N];
118 }
119 
120 void Diagnostic::define(const File &file, IO_Type default_type) const {
121  this->define_impl(file, default_type);
122 }
123 
124 //! Define NetCDF variables corresponding to a diagnostic quantity.
125 void Diagnostic::define_impl(const File &file, IO_Type default_type) const {
126  for (const auto &v : m_vars) {
127  io::define_spatial_variable(v, *m_grid, file, default_type);
128  }
129 }
130 
131 //! \brief A method for setting common variable attributes.
132 void Diagnostic::set_attrs(const std::string &long_name,
133  const std::string &standard_name,
134  const std::string &units,
135  const std::string &glaciological_units,
136  unsigned int N) {
137  if (N >= m_vars.size()) {
139  "N (%d) >= m_dof (%d)", N,
140  static_cast<int>(m_vars.size()));
141  }
142 
143  m_vars[N]["pism_intent"] = "diagnostic";
144 
145  m_vars[N]["long_name"] = long_name;
146 
147  m_vars[N]["standard_name"] = standard_name;
148 
149  if (units == glaciological_units) {
150  // No unit conversion needed to write data, so we assume that we don't need to check
151  // if units are valid. This is needed to be able to set units like "Pa s^(1/3)", which
152  // are correct (ice hardness with the Glen exponent n=3) but are not supported by
153  // UDUNITS.
154  //
155  // Also note that this automatically sets glaciological_units.
156  m_vars[N].set_units_without_validation(units);
157  } else {
158  m_vars[N]["units"] = units;
159 
160  if (not (m_config->get_flag("output.use_MKS") or glaciological_units.empty())) {
161  m_vars[N]["glaciological_units"] = glaciological_units;
162  }
163  }
164 }
165 
167  // use the name of the first variable
168  std::vector<std::string> names;
169  for (const auto &v : m_vars) {
170  names.push_back(v.get_name());
171  }
172  std::string all_names = join(names, ",");
173 
174  m_grid->ctx()->log()->message(3, "- Computing %s...\n", all_names.c_str());
175  IceModelVec::Ptr result = this->compute_impl();
176  m_grid->ctx()->log()->message(3, "- Done computing %s.\n", all_names.c_str());
177 
178  return result;
179 }
180 
181 TSDiagnostic::TSDiagnostic(IceGrid::ConstPtr grid, const std::string &name)
182  : m_grid(grid),
183  m_config(grid->ctx()->config()),
184  m_sys(grid->ctx()->unit_system()),
185  m_time_name(grid->ctx()->config()->get_string("time.dimension_name")),
186  m_variable(name, m_sys),
187  m_dimension(m_time_name, m_sys),
188  m_time_bounds(m_time_name + "_bounds", m_sys) {
189 
190  m_current_time = 0;
191  m_start = 0;
192 
193  m_buffer_size = static_cast<size_t>(m_config->get_number("output.timeseries.buffer_size"));
194 
195  m_variable["ancillary_variables"] = name + "_aux";
196 
197  m_dimension["calendar"] = m_grid->ctx()->time()->calendar();
198  m_dimension["long_name"] = "time";
199  m_dimension["axis"] = "T";
200  m_dimension["units"] = m_grid->ctx()->time()->units_string();
201 }
202 
204  flush();
205 }
206 
207 void TSDiagnostic::set_units(const std::string &units,
208  const std::string &glaciological_units) {
209  m_variable["units"] = units;
210 
211  if (not m_config->get_flag("output.use_MKS")) {
212  m_variable["glaciological_units"] = glaciological_units;
213  }
214 }
215 
217  : TSDiagnostic(grid, name) {
218  // empty
219 }
220 
222  : TSDiagnostic(grid, name), m_accumulator(0.0), m_v_previous(0.0), m_v_previous_set(false) {
223  // empty
224 }
225 
227  : TSRateDiagnostic(grid, name) {
228  // empty
229 }
230 
231 void TSSnapshotDiagnostic::evaluate(double t0, double t1, double v) {
232 
233  // skip times before the beginning of this time step
234  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] < t0) {
235  m_current_time += 1;
236  }
237 
238  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] <= t1) {
239  const unsigned int k = m_current_time;
240  m_current_time += 1;
241 
242  // skip the first time: it defines the beginning of a reporting interval
243  if (k == 0) {
244  continue;
245  }
246 
247  const double t_s = (*m_requested_times)[k - 1];
248  const double t_e = (*m_requested_times)[k];
249 
250  // store computed data in the buffer
251  {
252  m_time.push_back(t_e);
253  m_values.push_back(v);
254  m_bounds.push_back(t_s);
255  m_bounds.push_back(t_e);
256  }
257  }
258 }
259 
260 void TSRateDiagnostic::evaluate(double t0, double t1, double change) {
261  static const double epsilon = 1e-4; // seconds
262  assert(t1 > t0);
263 
264  // skip times before and including the beginning of this time step
265  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] <= t0) {
266  m_current_time += 1;
267  }
268 
269  // number of requested times in this time step
270  unsigned int N = 0;
271 
272  // loop through requested times that are within this time step
273  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] <= t1) {
274  const unsigned int k = m_current_time;
275  m_current_time += 1;
276 
277  N += 1;
278 
279  // skip the first time: it defines the beginning of a reporting interval
280  if (k == 0) {
281  continue;
282  }
283 
284  const double t_s = (*m_requested_times)[k - 1];
285  const double t_e = (*m_requested_times)[k];
286 
287  double rate = 0.0;
288  if (N == 1) {
289  // it is the right end-point of the first reporting interval in this time step: count the
290  // contribution from the last time step plus the one from the beginning of this time step
291  const double
292  total_change = m_accumulator + change * (t_e - t0) / (t1 - t0);
293  const double dt = t_e - t_s;
294 
295  rate = total_change / dt;
296 
297  } else {
298  // this reporting interval is completely contained within the time step, so the rate of change
299  // does not depend on its length
300  rate = change / (t1 - t0);
301  }
302 
303  // store computed data in the buffer
304  {
305  m_time.push_back(t_e);
306  m_values.push_back(rate);
307  m_bounds.push_back(t_s);
308  m_bounds.push_back(t_e);
309  }
310 
311  m_accumulator = 0.0;
312  }
313 
314  if (N == 0) {
315  // if this time step contained no requested times we need to add the whole change to the
316  // accumulator
317  m_accumulator += change;
318  } else {
319  // if this time step contained some requested times we need to add the change since the last one
320  // to the accumulator
321  const double dt = t1 - (*m_requested_times)[m_current_time - 1];
322  if (dt > epsilon) {
323  m_accumulator += change * (dt / (t1 - t0));
324  }
325  }
326 }
327 
328 void TSDiagnostic::update(double t0, double t1) {
329  this->update_impl(t0, t1);
330 }
331 
332 void TSSnapshotDiagnostic::update_impl(double t0, double t1) {
333  static const double epsilon = 1e-4; // seconds
334 
335  if (fabs(t1 - t0) < epsilon) {
336  return;
337  }
338 
339  assert(t1 > t0);
340 
341  evaluate(t0, t1, this->compute());
342 }
343 
344 void TSRateDiagnostic::update_impl(double t0, double t1) {
345  const double v = this->compute();
346 
347  if (m_v_previous_set) {
348  assert(t1 > t0);
349  evaluate(t0, t1, v - m_v_previous);
350  }
351 
352  m_v_previous = v;
353  m_v_previous_set = true;
354 }
355 
356 void TSFluxDiagnostic::update_impl(double t0, double t1) {
357  static const double epsilon = 1e-4; // seconds
358 
359  if (fabs(t1 - t0) < epsilon) {
360  return;
361  }
362 
363  assert(t1 > t0);
364 
365  evaluate(t0, t1, this->compute());
366 }
367 
368 void TSDiagnostic::define(const File &file) const {
369  auto time_name = m_config->get_string("time.dimension_name");
370  io::define_timeseries(m_variable, time_name, file, PISM_DOUBLE);
371  io::define_time_bounds(m_time_bounds, time_name, "nv", file, PISM_DOUBLE);
372 }
373 
375 
376  if (m_time.empty()) {
377  return;
378  }
379 
380  std::string dimension_name = m_dimension.get_name();
381 
382  File file(m_grid->com, m_output_filename, PISM_NETCDF3, PISM_READWRITE); // OK to use netcdf3
383 
384  unsigned int len = file.dimension_length(dimension_name);
385 
386  if (len > 0) {
387  double last_time = vector_max(file.read_dimension(dimension_name));
388 
389  if (last_time < m_time.front()) {
390  m_start = len;
391  }
392  }
393 
394  auto time_name = m_config->get_string("time.dimension_name");
395 
396  if (len == m_start) {
397  if (not file.find_variable(m_dimension.get_name())) {
398  io::define_timeseries(m_dimension, time_name, file, PISM_DOUBLE);
399  }
401 
402  if (not file.find_variable(m_time_bounds.get_name())) {
403  io::define_time_bounds(m_time_bounds, time_name, "nv", file, PISM_DOUBLE);
404  }
406  }
407 
408  if (not file.find_variable(m_variable.get_name())) {
409  io::define_timeseries(m_variable, time_name, file, PISM_DOUBLE);
410  }
412 
413  m_start += m_time.size();
414 
415  {
416  m_time.clear();
417  m_bounds.clear();
418  m_values.clear();
419  }
420 }
421 
422 void TSDiagnostic::init(const File &output_file,
423  std::shared_ptr<std::vector<double>> requested_times) {
424  m_output_filename = output_file.filename();
425 
426  m_requested_times = std::move(requested_times);
427 
428  // Get the number of records in the file (for appending):
429  m_start = output_file.dimension_length(m_dimension.get_name());
430 }
431 
433  return m_variable;
434 }
435 
436 } // end of namespace pism
void write_state(const File &output) const
Definition: Diagnostic.cc:89
virtual void write_state_impl(const File &output) const
Definition: Diagnostic.cc:104
IceModelVec::Ptr compute() const
Compute a diagnostic quantity and return a pointer to a newly-allocated IceModelVec.
Definition: Diagnostic.cc:166
virtual void reset_impl()
Definition: Diagnostic.cc:51
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:108
double to_internal(double x) const
Definition: Diagnostic.cc:58
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:112
virtual void define_state_impl(const File &output) const
Definition: Diagnostic.cc:99
void init(const File &input, unsigned int time)
Definition: Diagnostic.cc:81
virtual void update_impl(double dt)
Definition: Diagnostic.cc:42
void define(const File &file, IO_Type default_type) const
Definition: Diagnostic.cc:120
virtual IceModelVec::Ptr compute_impl() const =0
IceGrid::ConstPtr m_grid
the grid
Definition: Diagnostic.hh:106
double to_external(double x) const
Definition: Diagnostic.cc:69
void define_state(const File &output) const
Definition: Diagnostic.cc:85
virtual void define_impl(const File &file, IO_Type default_type) const
Define NetCDF variables corresponding to a diagnostic quantity.
Definition: Diagnostic.cc:125
unsigned int n_variables() const
Get the number of NetCDF variables corresponding to a diagnostic quantity.
Definition: Diagnostic.cc:77
void update(double dt)
Definition: Diagnostic.cc:38
Diagnostic(IceGrid::ConstPtr g)
Definition: Diagnostic.cc:30
virtual void init_impl(const File &input, unsigned int time)
Definition: Diagnostic.cc:93
void set_attrs(const std::string &long_name, const std::string &standard_name, const std::string &units, const std::string &glaciological_units, unsigned int N=0)
A method for setting common variable attributes.
Definition: Diagnostic.cc:132
SpatialVariableMetadata & metadata(unsigned int N=0)
Get a metadata object corresponding to variable number N.
Definition: Diagnostic.cc:110
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:110
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:364
std::string filename() const
Definition: File.cc:310
unsigned int dimension_length(const std::string &name) const
Get the length of a dimension.
Definition: File.cc:457
std::vector< double > read_dimension(const std::string &name) const
Get dimension data (a coordinate variable).
Definition: File.cc:596
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
std::shared_ptr< IceModelVec > Ptr
Definition: iceModelVec.hh:206
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
virtual ~TSDiagnostic()
Definition: Diagnostic.cc:203
std::vector< double > m_values
Definition: Diagnostic.hh:330
std::shared_ptr< std::vector< double > > m_requested_times
requested times
Definition: Diagnostic.hh:333
VariableMetadata m_dimension
Definition: Diagnostic.hh:324
unsigned int m_current_time
index into m_times
Definition: Diagnostic.hh:335
unsigned int m_start
starting index used when flushing the buffer
Definition: Diagnostic.hh:341
VariableMetadata m_time_bounds
Definition: Diagnostic.hh:325
std::vector< double > m_bounds
Definition: Diagnostic.hh:329
size_t m_buffer_size
size of the buffer used to store data
Definition: Diagnostic.hh:343
VariableMetadata m_variable
Definition: Diagnostic.hh:323
std::string m_output_filename
Definition: Diagnostic.hh:339
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:316
std::vector< double > m_time
Definition: Diagnostic.hh:328
IceGrid::ConstPtr m_grid
the grid
Definition: Diagnostic.hh:314
virtual double compute()=0
TSDiagnostic(IceGrid::ConstPtr g, const std::string &name)
Definition: Diagnostic.cc:181
void set_units(const std::string &units, const std::string &glaciological_units)
Definition: Diagnostic.cc:207
const VariableMetadata & metadata() const
Definition: Diagnostic.cc:432
virtual void update_impl(double t0, double t1)=0
void define(const File &file) const
Definition: Diagnostic.cc:368
void init(const File &output_file, std::shared_ptr< std::vector< double >> requested_times)
Definition: Diagnostic.cc:422
void update(double t0, double t1)
Definition: Diagnostic.cc:328
PISM's scalar time-series diagnostics.
Definition: Diagnostic.hh:278
void update_impl(double t0, double t1)
Definition: Diagnostic.cc:356
TSFluxDiagnostic(IceGrid::ConstPtr g, const std::string &name)
Definition: Diagnostic.cc:226
double m_accumulator
accumulator of changes (used to compute rates of change)
Definition: Diagnostic.hh:371
void update_impl(double t0, double t1)
Definition: Diagnostic.cc:344
void evaluate(double t0, double t1, double change)
Definition: Diagnostic.cc:260
double m_v_previous
last two values, used to compute the change during a time step
Definition: Diagnostic.hh:377
TSRateDiagnostic(IceGrid::ConstPtr g, const std::string &name)
Definition: Diagnostic.cc:221
Scalar diagnostic reporting the rate of change of a quantity modeled by PISM.
Definition: Diagnostic.hh:366
TSSnapshotDiagnostic(IceGrid::ConstPtr g, const std::string &name)
Definition: Diagnostic.cc:216
void evaluate(double t0, double t1, double v)
Definition: Diagnostic.cc:231
void update_impl(double t0, double t1)
Definition: Diagnostic.cc:332
std::string get_name() const
#define PISM_ERROR_LOCATION
void define_timeseries(const VariableMetadata &var, const std::string &dimension_name, const File &file, IO_Type nctype)
Define a NetCDF variable corresponding to a time-series.
Definition: io_helpers.cc:1026
void define_spatial_variable(const SpatialVariableMetadata &var, const IceGrid &grid, const File &file, IO_Type default_type)
Define a NetCDF variable corresponding to a VariableMetadata object.
Definition: io_helpers.cc:558
void write_time_bounds(const File &file, const VariableMetadata &metadata, size_t t_start, const std::vector< double > &data)
Definition: io_helpers.cc:1253
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:1147
void write_timeseries(const File &file, const VariableMetadata &metadata, size_t t_start, const std::vector< double > &data)
Write a time-series data to a file.
Definition: io_helpers.cc:1119
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:72
IO_Type
Definition: IO_Flags.hh:31
@ PISM_DOUBLE
Definition: IO_Flags.hh:38
static const double k
Definition: exactTestP.cc:45
@ PISM_NETCDF3
Definition: IO_Flags.hh:41
@ PISM_READWRITE
open an existing file for reading and writing
Definition: IO_Flags.hh:51
std::string join(const std::vector< std::string > &strings, const std::string &separator)
Concatenate strings, inserting separator between elements.
double vector_max(const std::vector< double > &input)