tDiagnostic.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tDiagnostic.hh (12493B)
---
1 // Copyright (C) 2010--2019 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 __Diagnostic_hh
20 #define __Diagnostic_hh
21
22 #include <memory>
23 #include <map>
24 #include <string>
25
26 #include "VariableMetadata.hh"
27 #include "Timeseries.hh" // inline code and a member of TSDiagnostic
28 #include "IceGrid.hh"
29 #include "ConfigInterface.hh"
30 #include "iceModelVec.hh"
31 #include "pism/util/error_handling.hh"
32 #include "pism/util/io/File.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:
61 Diagnostic(IceGrid::ConstPtr g);
62 virtual ~Diagnostic();
63
64 typedef std::shared_ptr<Diagnostic> Ptr;
65
66 static Ptr wrap(const IceModelVec2S &input);
67 static Ptr wrap(const IceModelVec2V &input);
68
69 void update(double dt);
70 void reset();
71
72 //! @brief Compute a diagnostic quantity and return a pointer to a newly-allocated IceModelVec.
73 IceModelVec::Ptr compute() const;
74
75 unsigned int n_variables() const;
76
77 SpatialVariableMetadata& metadata(unsigned int N = 0);
78
79 void define(const File &file, IO_Type default_type) const;
80
81 void init(const File &input, unsigned int time);
82 void define_state(const File &output) const;
83 void write_state(const File &output) const;
84 protected:
85 virtual void define_impl(const File &file, IO_Type default_type) const;
86 virtual void init_impl(const File &input, unsigned int time);
87 virtual void define_state_impl(const File &output) const;
88 virtual void write_state_impl(const File &output) const;
89
90 void set_attrs(const std::string &long_name,
91 const std::string &standard_name,
92 const std::string &units,
93 const std::string &glaciological_units,
94 unsigned int N = 0);
95
96 virtual void update_impl(double dt);
97 virtual void reset_impl();
98
99 virtual IceModelVec::Ptr compute_impl() const = 0;
100
101 double to_internal(double x) const;
102 double to_external(double x) const;
103
104 //! the grid
105 IceGrid::ConstPtr m_grid;
106 //! the unit system
107 const units::System::Ptr m_sys;
108 //! Configuration flags and parameters
109 const Config::ConstPtr m_config;
110 //! metadata corresponding to NetCDF variables
111 std::vector<SpatialVariableMetadata> m_vars;
112 //! fill value (used often enough to justify storing it)
113 double m_fill_value;
114 };
115
116 typedef std::map<std::string, Diagnostic::Ptr> DiagnosticList;
117
118 /*!
119 * Helper template wrapping quantities with dedicated storage in diagnostic classes.
120 *
121 * Note: Make sure that that created diagnostics don't outlast fields that they wrap (or you'll have
122 * dangling pointers).
123 */
124 template<class T>
125 class DiagWithDedicatedStorage : public Diagnostic {
126 public:
127 DiagWithDedicatedStorage(const T &input)
128 : Diagnostic(input.grid()),
129 m_input(input)
130 {
131 for (unsigned int j = 0; j < input.ndof(); ++j) {
132 m_vars.push_back(input.metadata(j));
133 }
134 }
135 protected:
136 IceModelVec::Ptr compute_impl() const {
137 typename T::Ptr result(new T(m_input.grid(), "unnamed", WITHOUT_GHOSTS));
138 result->set_name(m_input.get_name());
139 for (unsigned int k = 0; k < m_vars.size(); ++k) {
140 result->metadata(k) = m_vars[k];
141 }
142
143 result->copy_from(m_input);
144
145 return result;
146 }
147 const T &m_input;
148 };
149
150 //! A template derived from Diagnostic, adding a "Model".
151 template <class Model>
152 class Diag : public Diagnostic {
153 public:
154 Diag(const Model *m)
155 : Diagnostic(m->grid()), model(m) {}
156 protected:
157 const Model *model;
158 };
159
160 /*!
161 * Report a time-averaged rate of change of a quantity by accumulating changes over several time
162 * steps.
163 */
164 template<class M>
165 class DiagAverageRate : public Diag<M>
166 {
167 public:
168
169 enum InputKind {TOTAL_CHANGE = 0, RATE = 1};
170
171 DiagAverageRate(const M *m, const std::string &name, InputKind kind)
172 : Diag<M>(m),
173 m_factor(1.0),
174 m_input_kind(kind),
175 m_accumulator(Diagnostic::m_grid, name + "_accumulator", WITHOUT_GHOSTS),
176 m_interval_length(0.0),
177 m_time_since_reset(name + "_time_since_reset",
178 Diagnostic::m_config->get_string("time.dimension_name"),
179 Diagnostic::m_sys) {
180
181 m_time_since_reset.set_string("units", "seconds");
182 m_time_since_reset.set_string("long_name",
183 "time since " + m_accumulator.get_name() +
184 " was reset to 0");
185
186 m_accumulator.metadata().set_string("long_name",
187 "accumulator for the " + name + " diagnostic");
188
189 m_accumulator.set(0.0);
190 }
191 protected:
192 void init_impl(const File &input, unsigned int time) {
193 if (input.find_variable(m_accumulator.get_name())) {
194 m_accumulator.read(input, time);
195 } else {
196 m_accumulator.set(0.0);
197 }
198
199 if (input.find_variable(m_time_since_reset.get_name())) {
200 input.read_variable(m_time_since_reset.get_name(),
201 {time}, {1}, // start, count
202 &m_interval_length);
203 } else {
204 m_interval_length = 0.0;
205 }
206 }
207
208 void define_state_impl(const File &output) const {
209 m_accumulator.define(output);
210 io::define_timeseries(m_time_since_reset, output, PISM_DOUBLE);
211 }
212
213 void write_state_impl(const File &output) const {
214 m_accumulator.write(output);
215
216 const unsigned int
217 time_length = output.dimension_length(m_time_since_reset.get_dimension_name()),
218 t_start = time_length > 0 ? time_length - 1 : 0;
219 io::write_timeseries(output, m_time_since_reset, t_start, m_interval_length, PISM_DOUBLE);
220 }
221
222 virtual void update_impl(double dt) {
223 // Here the "factor" is used to convert units (from m to kg m-2, for example) and (possibly)
224 // integrate over the time integral using the rectangle method.
225
226 double factor = m_factor * (m_input_kind == TOTAL_CHANGE ? 1.0 : dt);
227
228 m_accumulator.add(factor, this->model_input());
229
230 m_interval_length += dt;
231 }
232
233 virtual void reset_impl() {
234 m_accumulator.set(0.0);
235 m_interval_length = 0.0;
236 }
237
238 virtual IceModelVec::Ptr compute_impl() const {
239 IceModelVec2S::Ptr result(new IceModelVec2S(Diagnostic::m_grid,
240 "diagnostic", WITHOUT_GHOSTS));
241 result->metadata(0) = Diagnostic::m_vars.at(0);
242
243 if (m_interval_length > 0.0) {
244 result->copy_from(m_accumulator);
245 result->scale(1.0 / m_interval_length);
246 } else {
247 result->set(Diagnostic::to_internal(Diagnostic::m_fill_value));
248 }
249
250 return result;
251 }
252 protected:
253 // constants initialized in the constructor
254 double m_factor;
255 InputKind m_input_kind;
256 // the state (read from and written to files)
257 IceModelVec2S m_accumulator;
258 // length of the reporting interval, accumulated along with the cumulative quantity
259 double m_interval_length;
260 TimeseriesMetadata m_time_since_reset;
261
262 // it should be enough to implement the constructor and this method
263 virtual const IceModelVec2S& model_input() {
264 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no default implementation");
265 }
266 };
267
268 //! @brief PISM's scalar time-series diagnostics.
269 class TSDiagnostic {
270 public:
271 typedef std::shared_ptr<TSDiagnostic> Ptr;
272
273 TSDiagnostic(IceGrid::ConstPtr g, const std::string &name);
274 virtual ~TSDiagnostic();
275
276 void update(double t0, double t1);
277
278 void flush();
279
280 void init(const File &output_file,
281 std::shared_ptr<std::vector<double>> requested_times);
282
283 const VariableMetadata &metadata() const;
284
285 void define(const File &file) const;
286
287 protected:
288 virtual void update_impl(double t0, double t1) = 0;
289
290 /*!
291 * Compute the diagnostic. Regular (snapshot) quantity should be computed here; for rates of
292 * change, compute() should return the total change during the time step from t0 to t1. The rate
293 * itself is computed in evaluate_rate().
294 */
295 virtual double compute() = 0;
296
297 /*!
298 * Set internal (MKS) and "glaciological" units.
299 *
300 * glaciological_units is ignored if output.use_MKS is set.
301 */
302 void set_units(const std::string &units, const std::string &glaciological_units);
303
304 //! the grid
305 IceGrid::ConstPtr m_grid;
306 //! Configuration flags and parameters
307 const Config::ConstPtr m_config;
308 //! the unit system
309 const units::System::Ptr m_sys;
310
311 //! time series object used to store computed values and metadata
312 Timeseries m_ts;
313
314 //! requested times
315 std::shared_ptr<std::vector<double>> m_times;
316 //! index into m_times
317 unsigned int m_current_time;
318
319 //! the name of the file to save to (stored here because it is used by flush(), which is called
320 //! from update())
321 std::string m_output_filename;
322 //! starting index used when flushing the buffer
323 unsigned int m_start;
324 //! size of the buffer used to store data
325 size_t m_buffer_size;
326 };
327
328 typedef std::map<std::string, TSDiagnostic::Ptr> TSDiagnosticList;
329
330 //! Scalar diagnostic reporting a snapshot of a quantity modeled by PISM.
331 /*!
332 * The method compute() should return the instantaneous "snapshot" value.
333 */
334 class TSSnapshotDiagnostic : public TSDiagnostic {
335 public:
336 TSSnapshotDiagnostic(IceGrid::ConstPtr g, const std::string &name);
337 private:
338 void update_impl(double t0, double t1);
339 void evaluate(double t0, double t1, double v);
340 };
341
342 //! Scalar diagnostic reporting the rate of change of a quantity modeled by PISM.
343 /*!
344 * The rate of change is averaged in time over reporting intervals.
345 *
346 * The method compute() should return the instantaneous "snapshot" value of a quantity.
347 */
348 class TSRateDiagnostic : public TSDiagnostic {
349 public:
350 TSRateDiagnostic(IceGrid::ConstPtr g, const std::string &name);
351 protected:
352 //! accumulator of changes (used to compute rates of change)
353 double m_accumulator;
354 void evaluate(double t0, double t1, double change);
355 private:
356 void update_impl(double t0, double t1);
357
358 //! last two values, used to compute the change during a time step
359 double m_v_previous;
360 bool m_v_previous_set;
361 };
362
363 //! Scalar diagnostic reporting a "flux".
364 /*!
365 * The flux is averaged over reporting intervals.
366 *
367 * The method compute() should return the change due to a flux over a time step.
368 *
369 * Fluxes can be computed using TSRateDiagnostic, but that would require keeping track of the total
370 * change due to a flux. It is possible for the magnitude of the total change to grow indefinitely,
371 * leading to the loss of precision; this is why we use changes over individual time steps instead.
372 *
373 * (The total change due to a flux can grow in magnitude even it the amount does not change. For
374 * example: if calving removes as much ice as we have added due to the SMB, the total mass is
375 * constant, but total SMB will grow.)
376 */
377 class TSFluxDiagnostic : public TSRateDiagnostic {
378 public:
379 TSFluxDiagnostic(IceGrid::ConstPtr g, const std::string &name);
380 private:
381 void update_impl(double t0, double t1);
382 };
383
384 template <class D, class M>
385 class TSDiag : public D {
386 public:
387 TSDiag(const M *m, const std::string &name)
388 : D(m->grid(), name), model(m) {
389 }
390 protected:
391 const M *model;
392 };
393
394 } // end of namespace pism
395
396 #endif /* __Diagnostic_hh */