tDiagnostic.cc - 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.cc (11153B)
---
1 /* Copyright (C) 2015, 2016, 2017, 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
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
27 namespace pism {
28
29 Diagnostic::Ptr Diagnostic::wrap(const IceModelVec2S &input) {
30 return Ptr(new DiagWithDedicatedStorage<IceModelVec2S>(input));
31 }
32
33 Diagnostic::Ptr Diagnostic::wrap(const IceModelVec2V &input) {
34 return Ptr(new DiagWithDedicatedStorage<IceModelVec2V>(input));
35 }
36
37 Diagnostic::Diagnostic(IceGrid::ConstPtr g)
38 : m_grid(g),
39 m_sys(g->ctx()->unit_system()),
40 m_config(g->ctx()->config()),
41 m_fill_value(m_config->get_number("output.fill_value")) {
42 // empty
43 }
44
45 Diagnostic::~Diagnostic() {
46 // empty
47 }
48
49 void Diagnostic::update(double dt) {
50 this->update_impl(dt);
51 }
52
53 void Diagnostic::update_impl(double dt) {
54 (void) dt;
55 // empty
56 }
57
58 void Diagnostic::reset() {
59 this->reset_impl();
60 }
61
62 void Diagnostic::reset_impl() {
63 // empty
64 }
65
66 /*!
67 * Convert from external (output) units to internal units.
68 */
69 double Diagnostic::to_internal(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, out, in);
74 }
75
76 /*!
77 * Convert from internal to external (output) units.
78 */
79 double Diagnostic::to_external(double x) const {
80 std::string
81 out = m_vars.at(0).get_string("glaciological_units"),
82 in = m_vars.at(0).get_string("units");
83 return convert(m_sys, x, in, out);
84 }
85
86 //! Get the number of NetCDF variables corresponding to a diagnostic quantity.
87 unsigned int Diagnostic::n_variables() const {
88 return m_vars.size();
89 }
90
91 void Diagnostic::init(const File &input, unsigned int time) {
92 this->init_impl(input, time);
93 }
94
95 void Diagnostic::define_state(const File &output) const {
96 this->define_state_impl(output);
97 }
98
99 void Diagnostic::write_state(const File &output) const {
100 this->write_state_impl(output);
101 }
102
103 void Diagnostic::init_impl(const File &input, unsigned int time) {
104 (void) input;
105 (void) time;
106 // empty
107 }
108
109 void Diagnostic::define_state_impl(const File &output) const {
110 (void) output;
111 // empty
112 }
113
114 void Diagnostic::write_state_impl(const File &output) const {
115 (void) output;
116 // empty
117 }
118
119 //! Get a metadata object corresponding to variable number N.
120 SpatialVariableMetadata& Diagnostic::metadata(unsigned int N) {
121 if (N >= m_vars.size()) {
122 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
123 "variable metadata index %d is out of bounds",
124 N);
125 }
126
127 return m_vars[N];
128 }
129
130 void Diagnostic::define(const File &file, IO_Type default_type) const {
131 this->define_impl(file, default_type);
132 }
133
134 //! Define NetCDF variables corresponding to a diagnostic quantity.
135 void Diagnostic::define_impl(const File &file, IO_Type default_type) const {
136 for (auto &v : m_vars) {
137 io::define_spatial_variable(v, *m_grid, file, default_type);
138 }
139 }
140
141 //! \brief A method for setting common variable attributes.
142 void Diagnostic::set_attrs(const std::string &long_name,
143 const std::string &standard_name,
144 const std::string &units,
145 const std::string &glaciological_units,
146 unsigned int N) {
147 if (N >= m_vars.size()) {
148 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
149 "N (%d) >= m_dof (%d)", N, (int)m_vars.size());
150 }
151
152 m_vars[N].set_string("pism_intent", "diagnostic");
153
154 m_vars[N].set_string("long_name", long_name);
155
156 m_vars[N].set_string("standard_name", standard_name);
157
158 m_vars[N].set_string("units", units);
159
160 if (not (m_config->get_flag("output.use_MKS") or glaciological_units.empty())) {
161 m_vars[N].set_string("glaciological_units", glaciological_units);
162 }
163 }
164
165 IceModelVec::Ptr Diagnostic::compute() const {
166 // use the name of the first variable
167 std::vector<std::string> names;
168 for (auto &v : m_vars) {
169 names.push_back(v.get_name());
170 }
171 std::string all_names = join(names, ",");
172
173 m_grid->ctx()->log()->message(3, "- Computing %s...\n", all_names.c_str());
174 IceModelVec::Ptr result = this->compute_impl();
175 m_grid->ctx()->log()->message(3, "- Done computing %s.\n", all_names.c_str());
176
177 return result;
178 }
179
180 TSDiagnostic::TSDiagnostic(IceGrid::ConstPtr g, const std::string &name)
181 : m_grid(g),
182 m_config(g->ctx()->config()),
183 m_sys(g->ctx()->unit_system()),
184 m_ts(*g, name, g->ctx()->config()->get_string("time.dimension_name")) {
185
186 m_current_time = 0;
187 m_start = 0;
188
189 m_buffer_size = (size_t)m_config->get_number("output.timeseries.buffer_size");
190
191 m_ts.variable().set_string("ancillary_variables", name + "_aux");
192
193 m_ts.dimension().set_string("calendar", m_grid->ctx()->time()->calendar());
194 m_ts.dimension().set_string("long_name", m_config->get_string("time.dimension_name"));
195 m_ts.dimension().set_string("axis", "T");
196 m_ts.dimension().set_string("units", m_grid->ctx()->time()->CF_units_string());
197 }
198
199 TSDiagnostic::~TSDiagnostic() {
200 flush();
201 }
202
203 void TSDiagnostic::set_units(const std::string &units,
204 const std::string &glaciological_units) {
205 m_ts.variable().set_string("units", units);
206
207 if (not m_config->get_flag("output.use_MKS")) {
208 m_ts.variable().set_string("glaciological_units", glaciological_units);
209 }
210 }
211
212 TSSnapshotDiagnostic::TSSnapshotDiagnostic(IceGrid::ConstPtr g, const std::string &name)
213 : TSDiagnostic(g, name) {
214 // empty
215 }
216
217 TSRateDiagnostic::TSRateDiagnostic(IceGrid::ConstPtr g, const std::string &name)
218 : TSDiagnostic(g, name), m_accumulator(0.0), m_v_previous(0.0), m_v_previous_set(false) {
219 // empty
220 }
221
222 TSFluxDiagnostic::TSFluxDiagnostic(IceGrid::ConstPtr g, const std::string &name)
223 : TSRateDiagnostic(g, name) {
224 // empty
225 }
226
227 void TSSnapshotDiagnostic::evaluate(double t0, double t1, double v) {
228
229 // skip times before the beginning of this time step
230 while (m_current_time < m_times->size() and (*m_times)[m_current_time] < t0) {
231 m_current_time += 1;
232 }
233
234 while (m_current_time < m_times->size() and (*m_times)[m_current_time] <= t1) {
235 const unsigned int k = m_current_time;
236 m_current_time += 1;
237
238 // skip the first time: it defines the beginning of a reporting interval
239 if (k == 0) {
240 continue;
241 }
242
243 const double
244 t_s = (*m_times)[k - 1],
245 t_e = (*m_times)[k];
246
247 m_ts.append(v, t_s, t_e);
248 }
249 }
250
251 void TSRateDiagnostic::evaluate(double t0, double t1, double change) {
252 static const double epsilon = 1e-4; // seconds
253 assert(t1 > t0);
254
255 // skip times before and including the beginning of this time step
256 while (m_current_time < m_times->size() and (*m_times)[m_current_time] <= t0) {
257 m_current_time += 1;
258 }
259
260 // number of requested times in this time step
261 unsigned int N = 0;
262
263 // loop through requested times that are within this time step
264 while (m_current_time < m_times->size() and (*m_times)[m_current_time] <= t1) {
265 const unsigned int k = m_current_time;
266 m_current_time += 1;
267
268 N += 1;
269
270 // skip the first time: it defines the beginning of a reporting interval
271 if (k == 0) {
272 continue;
273 }
274
275 const double
276 t_s = (*m_times)[k - 1],
277 t_e = (*m_times)[k];
278
279 double rate = 0.0;
280 if (N == 1) {
281 // it is the right end-point of the first reporting interval in this time step: count the
282 // contribution from the last time step plus the one from the beginning of this time step
283 const double
284 total_change = m_accumulator + change * (t_e - t0) / (t1 - t0),
285 dt = t_e - t_s;
286
287 rate = total_change / dt;
288
289 } else {
290 // this reporting interval is completely contained within the time step, so the rate of change
291 // does not depend on its length
292 rate = change / (t1 - t0);
293 }
294
295 m_ts.append(rate, t_s, t_e);
296 m_accumulator = 0.0;
297 }
298
299 if (N == 0) {
300 // if this time step contained no requested times we need to add the whole change to the
301 // accumulator
302 m_accumulator += change;
303 } else {
304 // if this time step contained some requested times we need to add the change since the last one
305 // to the accumulator
306 const double dt = t1 - (*m_times)[m_current_time - 1];
307 if (dt > epsilon) {
308 m_accumulator += change * (dt / (t1 - t0));
309 }
310 }
311 }
312
313 void TSDiagnostic::update(double t0, double t1) {
314 this->update_impl(t0, t1);
315 }
316
317 void TSSnapshotDiagnostic::update_impl(double t0, double t1) {
318 if (fabs(t1 - t0) < 1e-2) {
319 return;
320 }
321 assert(t1 > t0);
322 evaluate(t0, t1, this->compute());
323 }
324
325 void TSRateDiagnostic::update_impl(double t0, double t1) {
326 const double v = this->compute();
327
328 if (m_v_previous_set) {
329 assert(t1 > t0);
330 evaluate(t0, t1, v - m_v_previous);
331 }
332
333 m_v_previous = v;
334 m_v_previous_set = true;
335 }
336
337 void TSFluxDiagnostic::update_impl(double t0, double t1) {
338 if (fabs(t1 - t0) < 1e-2) {
339 return;
340 }
341 assert(t1 > t0);
342 evaluate(t0, t1, this->compute());
343 }
344
345 void TSDiagnostic::define(const File &file) const {
346 io::define_timeseries(m_ts.variable(), file, PISM_DOUBLE);
347 io::define_time_bounds(m_ts.bounds(), file, PISM_DOUBLE);
348 }
349
350 void TSDiagnostic::flush() {
351
352 if (m_ts.times().empty()) {
353 return;
354 }
355
356 std::string dimension_name = m_ts.dimension().get_name();
357
358 File file(m_grid->com, m_output_filename, PISM_NETCDF3, PISM_READWRITE); // OK to use netcdf3
359
360 unsigned int len = file.dimension_length(dimension_name);
361
362 if (len > 0) {
363 double last_time = vector_max(file.read_dimension(dimension_name));
364
365 if (last_time < m_ts.times().front()) {
366 m_start = len;
367 }
368 }
369
370 if (len == m_start) {
371 io::write_timeseries(file, m_ts.dimension(), m_start, m_ts.times());
372 io::write_time_bounds(file, m_ts.bounds(), m_start, m_ts.time_bounds());
373 }
374 io::write_timeseries(file, m_ts.variable(), m_start, m_ts.values());
375
376 m_start += m_ts.times().size();
377
378 m_ts.reset();
379 }
380
381 void TSDiagnostic::init(const File &output_file,
382 std::shared_ptr<std::vector<double>> requested_times) {
383 m_output_filename = output_file.filename();
384
385 m_times = requested_times;
386
387 // Get the number of records in the file (for appending):
388 m_start = output_file.dimension_length(m_ts.dimension().get_name());
389 }
390
391 const VariableMetadata &TSDiagnostic::metadata() const {
392 return m_ts.variable();
393 }
394
395 } // end of namespace pism