PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Diagnostic.hh
Go to the documentation of this file.
1 // Copyright (C) 2010--2021 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 #ifndef PISM_DIAGNOSTIC_HH
20 #define PISM_DIAGNOSTIC_HH
21 
22 #include <memory>
23 #include <map>
24 #include <string>
25 
26 #include "VariableMetadata.hh"
27 #include "IceGrid.hh"
28 #include "ConfigInterface.hh"
29 #include "iceModelVec.hh"
30 #include "pism/util/error_handling.hh"
31 #include "pism/util/io/File.hh"
32 #include "pism/util/IceModelVec2V.hh"
33 #include "pism/util/io/io_helpers.hh"
34 
35 namespace pism {
36 
37 //! @brief Class representing diagnostic computations in PISM.
38 /*!
39  * The main goal of this abstraction is to allow accessing metadata
40  * corresponding to a diagnostic quantity *before* it is computed.
41  *
42  * Another goal is to create an interface for computing diagnostics *without*
43  * knowing which PISM module is responsible for the computation.
44  *
45  * Technical note: to compute some diagnostic quantities we need access to
46  * protected members of classes. C++ forbids obtaining pointers to non-static
47  * methods of a class, but it is possible to define a (friend) function
48  *
49  * @code
50  * IceModelVec::Ptr compute_bar(Foo* model, ...);
51  * @endcode
52  *
53  * which is the same as creating a method `Foo::compute_bar()`, but you *can*
54  * get a pointer to it.
55  *
56  * Diagnostic creates a common interface for all these compute_bar
57  * functions.
58  */
59 class Diagnostic {
60 public:
62  virtual ~Diagnostic() = default;
63 
64  typedef std::shared_ptr<Diagnostic> Ptr;
65 
66  // defined below
67  template<typename T>
68  static Ptr wrap(const T &input);
69 
70  void update(double dt);
71  void reset();
72 
73  //! @brief Compute a diagnostic quantity and return a pointer to a newly-allocated IceModelVec.
74  IceModelVec::Ptr compute() const;
75 
76  unsigned int n_variables() const;
77 
78  SpatialVariableMetadata& metadata(unsigned int N = 0);
79 
80  void define(const File &file, IO_Type default_type) const;
81 
82  void init(const File &input, unsigned int time);
83  void define_state(const File &output) const;
84  void write_state(const File &output) const;
85 protected:
86  virtual void define_impl(const File &file, IO_Type default_type) const;
87  virtual void init_impl(const File &input, unsigned int time);
88  virtual void define_state_impl(const File &output) const;
89  virtual void write_state_impl(const File &output) const;
90 
91  void set_attrs(const std::string &long_name,
92  const std::string &standard_name,
93  const std::string &units,
94  const std::string &glaciological_units,
95  unsigned int N = 0);
96 
97  virtual void update_impl(double dt);
98  virtual void reset_impl();
99 
100  virtual IceModelVec::Ptr compute_impl() const = 0;
101 
102  double to_internal(double x) const;
103  double to_external(double x) const;
104 
105  //! the grid
107  //! the unit system
109  //! Configuration flags and parameters
111  //! metadata corresponding to NetCDF variables
112  std::vector<SpatialVariableMetadata> m_vars;
113  //! fill value (used often enough to justify storing it)
114  double m_fill_value;
115 };
116 
117 typedef std::map<std::string, Diagnostic::Ptr> DiagnosticList;
118 
119 /*!
120  * Helper template wrapping quantities with dedicated storage in diagnostic classes.
121  *
122  * Note: Make sure that that created diagnostics don't outlast fields that they wrap (or you'll have
123  * dangling pointers).
124  */
125 template<class T>
127 public:
128  DiagWithDedicatedStorage(const T &input)
129  : Diagnostic(input.grid()),
130  m_input(input)
131  {
132  for (unsigned int j = 0; j < input.ndof(); ++j) {
133  m_vars.emplace_back(input.metadata(j));
134  }
135  }
136 protected:
137 
139  auto result = duplicate(m_input);
140 
141  result->set_name(m_input.get_name());
142  for (unsigned int k = 0; k < m_vars.size(); ++k) {
143  result->metadata(k) = m_vars[k];
144  }
145 
146  result->copy_from(m_input);
147 
148  return result;
149  }
150 
151  const T &m_input;
152 };
153 
154 template<typename T>
156  return Ptr(new DiagWithDedicatedStorage<T>(input));
157 }
158 
159 //! A template derived from Diagnostic, adding a "Model".
160 template <class Model>
161 class Diag : public Diagnostic {
162 public:
163  Diag(const Model *m)
164  : Diagnostic(m->grid()), model(m) {}
165 protected:
166  const Model *model;
167 };
168 
169 /*!
170  * Report a time-averaged rate of change of a quantity by accumulating changes over several time
171  * steps.
172  */
173 template<class M>
174 class DiagAverageRate : public Diag<M>
175 {
176 public:
177 
178  enum InputKind {TOTAL_CHANGE = 0, RATE = 1};
179 
180  DiagAverageRate(const M *m, const std::string &name, InputKind kind)
181  : Diag<M>(m),
182  m_factor(1.0),
183  m_input_kind(kind),
184  m_accumulator(Diagnostic::m_grid, name + "_accumulator", WITHOUT_GHOSTS),
185  m_interval_length(0.0),
186  m_time_since_reset(name + "_time_since_reset", Diagnostic::m_sys) {
187 
188  m_time_since_reset["units"] = "seconds";
189  m_time_since_reset["long_name"] =
190  "time since " + m_accumulator.get_name() + " was reset to 0";
191 
192  m_accumulator.metadata()["long_name"] =
193  "accumulator for the " + name + " diagnostic";
194 
195  m_accumulator.set(0.0);
196  }
197 protected:
198  void init_impl(const File &input, unsigned int time) {
199  if (input.find_variable(m_accumulator.get_name())) {
200  m_accumulator.read(input, time);
201  } else {
202  m_accumulator.set(0.0);
203  }
204 
207  {time}, {1}, // start, count
209  } else {
210  m_interval_length = 0.0;
211  }
212  }
213 
214  void define_state_impl(const File &output) const {
215  m_accumulator.define(output);
217  Diagnostic::m_config->get_string("time.dimension_name"),
218  output, PISM_DOUBLE);
219  }
220 
221  void write_state_impl(const File &output) const {
222  m_accumulator.write(output);
223 
224  auto time_name = Diagnostic::m_config->get_string("time.dimension_name");
225 
226  unsigned int time_length = output.dimension_length(time_name);
227  unsigned int t_start = time_length > 0 ? time_length - 1 : 0;
229  }
230 
231  virtual void update_impl(double dt) {
232  // Here the "factor" is used to convert units (from m to kg m-2, for example) and (possibly)
233  // integrate over the time integral using the rectangle method.
234 
235  double factor = m_factor * (m_input_kind == TOTAL_CHANGE ? 1.0 : dt);
236 
237  m_accumulator.add(factor, this->model_input());
238 
239  m_interval_length += dt;
240  }
241 
242  virtual void reset_impl() {
243  m_accumulator.set(0.0);
244  m_interval_length = 0.0;
245  }
246 
247  virtual IceModelVec::Ptr compute_impl() const {
249  "diagnostic", WITHOUT_GHOSTS));
250  result->metadata(0) = Diagnostic::m_vars.at(0);
251 
252  if (m_interval_length > 0.0) {
253  result->copy_from(m_accumulator);
254  result->scale(1.0 / m_interval_length);
255  } else {
257  }
258 
259  return result;
260  }
261 
262  // constants initialized in the constructor
263  double m_factor;
265  // the state (read from and written to files)
267  // length of the reporting interval, accumulated along with the cumulative quantity
270 
271  // it should be enough to implement the constructor and this method
272  virtual const IceModelVec2S& model_input() {
273  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no default implementation");
274  }
275 };
276 
277 //! @brief PISM's scalar time-series diagnostics.
279 public:
280  typedef std::shared_ptr<TSDiagnostic> Ptr;
281 
282  TSDiagnostic(IceGrid::ConstPtr g, const std::string &name);
283  virtual ~TSDiagnostic();
284 
285  void update(double t0, double t1);
286 
287  void flush();
288 
289  void init(const File &output_file,
290  std::shared_ptr<std::vector<double>> requested_times);
291 
292  const VariableMetadata &metadata() const;
293 
294  void define(const File &file) const;
295 
296 protected:
297  virtual void update_impl(double t0, double t1) = 0;
298 
299  /*!
300  * Compute the diagnostic. Regular (snapshot) quantity should be computed here; for rates of
301  * change, compute() should return the total change during the time step from t0 to t1. The rate
302  * itself is computed in evaluate_rate().
303  */
304  virtual double compute() = 0;
305 
306  /*!
307  * Set internal (MKS) and "glaciological" units.
308  *
309  * glaciological_units is ignored if output.use_MKS is set.
310  */
311  void set_units(const std::string &units, const std::string &glaciological_units);
312 
313  //! the grid
315  //! Configuration flags and parameters
317  //! the unit system
319 
320  //! time series object used to store computed values and metadata
321  std::string m_time_name;
322 
326 
327  // buffer for diagnostic time series
328  std::vector<double> m_time;
329  std::vector<double> m_bounds;
330  std::vector<double> m_values;
331 
332  //! requested times
333  std::shared_ptr<std::vector<double>> m_requested_times;
334  //! index into m_times
335  unsigned int m_current_time;
336 
337  //! the name of the file to save to (stored here because it is used by flush(), which is called
338  //! from update())
339  std::string m_output_filename;
340  //! starting index used when flushing the buffer
341  unsigned int m_start;
342  //! size of the buffer used to store data
344 };
345 
346 typedef std::map<std::string, TSDiagnostic::Ptr> TSDiagnosticList;
347 
348 //! Scalar diagnostic reporting a snapshot of a quantity modeled by PISM.
349 /*!
350  * The method compute() should return the instantaneous "snapshot" value.
351  */
353 public:
354  TSSnapshotDiagnostic(IceGrid::ConstPtr g, const std::string &name);
355 private:
356  void update_impl(double t0, double t1);
357  void evaluate(double t0, double t1, double v);
358 };
359 
360 //! Scalar diagnostic reporting the rate of change of a quantity modeled by PISM.
361 /*!
362  * The rate of change is averaged in time over reporting intervals.
363  *
364  * The method compute() should return the instantaneous "snapshot" value of a quantity.
365  */
367 public:
368  TSRateDiagnostic(IceGrid::ConstPtr g, const std::string &name);
369 protected:
370  //! accumulator of changes (used to compute rates of change)
372  void evaluate(double t0, double t1, double change);
373 private:
374  void update_impl(double t0, double t1);
375 
376  //! last two values, used to compute the change during a time step
377  double m_v_previous;
379 };
380 
381 //! Scalar diagnostic reporting a "flux".
382 /*!
383  * The flux is averaged over reporting intervals.
384  *
385  * The method compute() should return the change due to a flux over a time step.
386  *
387  * Fluxes can be computed using TSRateDiagnostic, but that would require keeping track of the total
388  * change due to a flux. It is possible for the magnitude of the total change to grow indefinitely,
389  * leading to the loss of precision; this is why we use changes over individual time steps instead.
390  *
391  * (The total change due to a flux can grow in magnitude even it the amount does not change. For
392  * example: if calving removes as much ice as we have added due to the SMB, the total mass is
393  * constant, but total SMB will grow.)
394  */
396 public:
397  TSFluxDiagnostic(IceGrid::ConstPtr g, const std::string &name);
398 private:
399  void update_impl(double t0, double t1);
400 };
401 
402 template <class D, class M>
403 class TSDiag : public D {
404 public:
405  TSDiag(const M *m, const std::string &name)
406  : D(m->grid(), name), model(m) {
407  }
408 protected:
409  const M *model;
410 };
411 
412 } // end of namespace pism
413 
414 #endif /* PISM_DIAGNOSTIC_HH */
std::shared_ptr< const Config > ConstPtr
virtual void update_impl(double dt)
Definition: Diagnostic.hh:231
VariableMetadata m_time_since_reset
Definition: Diagnostic.hh:269
virtual void reset_impl()
Definition: Diagnostic.hh:242
void init_impl(const File &input, unsigned int time)
Definition: Diagnostic.hh:198
IceModelVec2S m_accumulator
Definition: Diagnostic.hh:266
virtual const IceModelVec2S & model_input()
Definition: Diagnostic.hh:272
virtual IceModelVec::Ptr compute_impl() const
Definition: Diagnostic.hh:247
DiagAverageRate(const M *m, const std::string &name, InputKind kind)
Definition: Diagnostic.hh:180
void write_state_impl(const File &output) const
Definition: Diagnostic.hh:221
void define_state_impl(const File &output) const
Definition: Diagnostic.hh:214
DiagWithDedicatedStorage(const T &input)
Definition: Diagnostic.hh:128
IceModelVec::Ptr compute_impl() const
Definition: Diagnostic.hh:138
const Model * model
Definition: Diagnostic.hh:166
Diag(const Model *m)
Definition: Diagnostic.hh:163
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:161
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
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
virtual void reset_impl()
Definition: Diagnostic.cc:51
virtual ~Diagnostic()=default
double m_fill_value
fill value (used often enough to justify storing it)
Definition: Diagnostic.hh:114
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
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:64
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
Class representing diagnostic computations in PISM.
Definition: Diagnostic.hh:59
void read_variable(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, double *ip) const
Definition: File.cc:740
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
unsigned int dimension_length(const std::string &name) const
Get the length of a dimension.
Definition: File.cc:457
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
void add(double alpha, const IceModelVec2S &x)
std::shared_ptr< IceModelVec2S > Ptr
Definition: iceModelVec.hh:341
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
std::shared_ptr< IceModelVec > Ptr
Definition: iceModelVec.hh:206
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void read(const std::string &filename, unsigned int time)
Definition: iceModelVec.cc:833
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
Definition: iceModelVec.cc:523
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
const std::string & get_name() const
Get the name of an IceModelVec object.
Definition: iceModelVec.cc:386
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).
TSDiag(const M *m, const std::string &name)
Definition: Diagnostic.hh:405
const M * model
Definition: Diagnostic.hh:409
virtual ~TSDiagnostic()
Definition: Diagnostic.cc:203
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:318
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::string m_time_name
time series object used to store computed values and metadata
Definition: Diagnostic.hh:321
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
std::shared_ptr< TSDiagnostic > Ptr
Definition: Diagnostic.hh:280
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
Scalar diagnostic reporting a "flux".
Definition: Diagnostic.hh:395
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
Scalar diagnostic reporting a snapshot of a quantity modeled by PISM.
Definition: Diagnostic.hh:352
std::string get_name() const
std::shared_ptr< System > Ptr
Definition: Units.hh:47
#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 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
IO_Type
Definition: IO_Flags.hh:31
@ PISM_DOUBLE
Definition: IO_Flags.hh:38
static const double g
Definition: exactTestP.cc:39
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:346
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
static const double k
Definition: exactTestP.cc:45
std::shared_ptr< IceModelVec2S > duplicate(const IceModelVec2S &source)
Definition: iceModelVec2.cc:55
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49