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