tTimeseries.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
---
tTimeseries.cc (9814B)
---
1 // Copyright (C) 2009--2020 Constantine Khroulev
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 #include <algorithm> // min_element and max_element
20 #include <gsl/gsl_interp.h>
21
22 #include "Timeseries.hh"
23 #include "pism_utilities.hh"
24 #include "IceGrid.hh"
25 #include "pism/util/io/File.hh"
26 #include "Time.hh"
27 #include "ConfigInterface.hh"
28
29 #include "error_handling.hh"
30 #include "io/io_helpers.hh"
31 #include "pism/util/Logger.hh"
32
33 namespace pism {
34
35 Timeseries::Timeseries(const IceGrid &g, const std::string &name, const std::string &dimension_name)
36 : m_dimension(dimension_name, dimension_name, g.ctx()->unit_system()),
37 m_variable(name, dimension_name, g.ctx()->unit_system()),
38 m_bounds(dimension_name + "_bounds", dimension_name, g.ctx()->unit_system())
39 {
40 private_constructor(g.com, dimension_name);
41 }
42
43 Timeseries::Timeseries(MPI_Comm c, units::System::Ptr unit_system,
44 const std::string & name, const std::string & dimension_name)
45 : m_dimension(dimension_name, dimension_name, unit_system),
46 m_variable(name, dimension_name, unit_system),
47 m_bounds(dimension_name + "_bounds", dimension_name, unit_system)
48 {
49 private_constructor(c, dimension_name);
50 }
51
52 void Timeseries::private_constructor(MPI_Comm c, const std::string &dimension_name) {
53 m_com = c;
54 m_dimension.set_string("bounds", dimension_name + "_bounds");
55
56 m_use_bounds = true;
57 }
58
59 //! Ensure that time bounds have the same units as the dimension.
60 void Timeseries::set_bounds_units() {
61 m_bounds.set_string("units", m_dimension.get_string("units"));
62 m_bounds.set_string("glaciological_units", m_dimension.get_string("glaciological_units"));
63 }
64
65 bool Timeseries::get_use_bounds() const {
66 return m_use_bounds;
67 }
68
69 void Timeseries::set_use_bounds(bool flag) {
70 m_use_bounds = flag;
71 }
72
73
74 //! Read timeseries data from a NetCDF file `filename`.
75 void Timeseries::read(const File &nc, const Time &time_manager, const Logger &log) {
76
77 std::string standard_name = m_variable.get_string("standard_name");
78
79 auto var = nc.find_variable(name(), standard_name);
80
81 if (not var.exists) {
82 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Can't find '%s' ('%s') in '%s'.\n",
83 name().c_str(), standard_name.c_str(),
84 nc.filename().c_str());
85 }
86
87 auto dims = nc.dimensions(var.name);
88
89 if (dims.size() != 1) {
90 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
91 "Variable '%s' in '%s' depends on %d dimensions,\n"
92 "but a time-series variable can only depend on 1 dimension.",
93 name().c_str(),
94 nc.filename().c_str(),
95 (int)dims.size());
96 }
97
98 auto time_name = dims[0];
99
100 TimeseriesMetadata tmp_dim = m_dimension;
101 tmp_dim.set_name(time_name);
102
103 io::read_timeseries(nc, tmp_dim, time_manager, log, m_time);
104
105 if (not is_increasing(m_time)) {
106 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
107 "dimension '%s' has to be strictly increasing (read from '%s').",
108 tmp_dim.get_name().c_str(), nc.filename().c_str());
109 }
110
111 std::string time_bounds_name = nc.read_text_attribute(time_name, "bounds");
112
113 if (not time_bounds_name.empty()) {
114 m_use_bounds = true;
115
116 set_bounds_units();
117 TimeBoundsMetadata tmp_bounds = m_bounds;
118 tmp_bounds.set_name(time_bounds_name);
119
120 tmp_bounds.set_string("units", tmp_dim.get_string("units"));
121
122 io::read_time_bounds(nc, tmp_bounds, time_manager, log, m_time_bounds);
123
124 // Time bounds override the time dimension read from the input file.
125 //
126 // NB: we use the right end point of each interval to be consistent with
127 // Timeseries::append()
128 for (unsigned int k = 0; k < m_time.size(); ++k) {
129 m_time[k] = m_time_bounds[2 * k + 1];
130 }
131 } else {
132 m_use_bounds = false;
133 }
134
135 io::read_timeseries(nc, m_variable, time_manager, log, m_values);
136
137 if (m_time.size() != m_values.size()) {
138 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "variables %s and %s in %s have different numbers of values.",
139 m_dimension.get_name().c_str(),
140 m_variable.get_name().c_str(),
141 nc.filename().c_str());
142 }
143
144 report_range(log);
145 }
146
147 //! \brief Report the range of a time-series stored in `values`.
148 void Timeseries::report_range(const Logger &log) {
149 double min, max;
150
151 // min_element and max_element return iterators; "*" is used to get
152 // the value corresponding to this iterator
153 min = *std::min_element(m_values.begin(), m_values.end());
154 max = *std::max_element(m_values.begin(), m_values.end());
155
156 units::Converter c(m_variable.unit_system(),
157 m_variable.get_string("units"),
158 m_variable.get_string("glaciological_units"));
159 min = c(min);
160 max = c(max);
161
162 std::string spacer(m_variable.get_name().size(), ' ');
163
164 log.message(2,
165 " FOUND %s / %-60s\n"
166 " %s \\ min,max = %9.3f,%9.3f (%s)\n",
167 m_variable.get_name().c_str(),
168 m_variable.get_string("long_name").c_str(), spacer.c_str(), min, max,
169 m_variable.get_string("glaciological_units").c_str());
170 }
171
172 //! Write timeseries data to a NetCDF file `filename`.
173 void Timeseries::write(const File &nc) const {
174 // write the dimensional variable; this call should go first
175 io::write_timeseries(nc, m_dimension, 0, m_time);
176 io::write_timeseries(nc, m_variable, 0, m_values);
177
178 if (m_use_bounds) {
179 io::write_time_bounds(nc, m_bounds, 0, m_time_bounds);
180 }
181 }
182
183 //! Clear storage.
184 void Timeseries::reset() {
185 m_time.clear();
186 m_values.clear();
187 m_time_bounds.clear();
188 }
189
190 /** Scale all values stored in this instance by `scaling_factor`.
191 *
192 * This is used to convert mass balance offsets from [m s-1] to [kg m-2 s-1].
193 *
194 * @param[in] scaling_factor multiplicative scaling factor
195 */
196 void Timeseries::scale(double scaling_factor) {
197 for (unsigned int i = 0; i < m_values.size(); ++i) {
198 m_values[i] *= scaling_factor;
199 }
200 }
201
202 //! Get a value of timeseries at time `t`.
203 /*! Returns the first value or the last value if t is out of range on the left
204 and right, respectively.
205
206 Uses time bounds if present (interpreting data as piecewise-constant) and
207 uses linear interpolation otherwise.
208 */
209 double Timeseries::operator()(double t) const {
210
211 if (m_use_bounds) {
212 // piecewise-constant case
213
214 size_t k = 0;
215 if (t < m_time[0]) {
216 k = 0;
217 } else if (t >= m_time.back()) {
218 k = m_time.size() - 1;
219 } else {
220 k = gsl_interp_bsearch(m_time.data(), t, 0, m_time.size()) + 1;
221 }
222
223 return m_values[k];
224 } else {
225 // piecewise-linear case
226
227 // extrapolation on the left
228 if (t < m_time[0]) {
229 return m_values[0];
230 }
231
232 size_t k = gsl_interp_bsearch(m_time.data(), t, 0, m_time.size());
233
234 // extrapolation on the right
235 if (k + 1 >= m_time.size()) {
236 return m_values.back();
237 }
238
239 double alpha = (t - m_time[k]) / (m_time[k + 1] - m_time[k]);
240
241 return m_values[k] * (1.0 - alpha) + m_values[k + 1] * alpha;
242 }
243 }
244
245 //! Get a value of timeseries by index.
246 /*!
247 Stops if the index is out of range.
248 */
249 double Timeseries::operator[](unsigned int j) const {
250
251 #if (Pism_DEBUG==1)
252 if (j >= m_values.size()) {
253 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Timeseries %s: operator[]: invalid argument: size=%d, index=%d",
254 m_variable.get_name().c_str(), (int)m_values.size(), j);
255 }
256 #endif
257
258 return m_values[j];
259 }
260
261 //! \brief Compute an average of a time-series over interval (t,t+dt) using
262 //! trapezoidal rule with N sub-intervals.
263 double Timeseries::average(double t, double dt, unsigned int N) const {
264 std::vector<double> V(N+1);
265
266 for (unsigned int i = 0; i < N+1; ++i) {
267 double t_i = t + (dt / N) * i;
268 V[i] = (*this)(t_i);
269 }
270
271 double sum = 0;
272 for (unsigned int i = 0; i < N; ++i) {
273 sum += V[i] + V[i+1];
274 }
275
276 return sum / (2.0*N);
277 }
278
279 //! Append a pair (t,v) to the timeseries.
280 void Timeseries::append(double v, double t0, double t1) {
281 m_time.push_back(t1);
282 m_values.push_back(v);
283 m_time_bounds.push_back(t0);
284 m_time_bounds.push_back(t1);
285 }
286
287 std::string Timeseries::name() const {
288 return m_variable.get_name();
289 }
290
291 TimeseriesMetadata& Timeseries::variable() {
292 return m_variable;
293 }
294
295 TimeseriesMetadata& Timeseries::dimension() {
296 return m_dimension;
297 }
298
299 TimeBoundsMetadata& Timeseries::bounds() {
300 return m_bounds;
301 }
302
303 const TimeseriesMetadata& Timeseries::variable() const {
304 return m_variable;
305 }
306
307 const TimeseriesMetadata& Timeseries::dimension() const {
308 return m_dimension;
309 }
310
311 const TimeBoundsMetadata& Timeseries::bounds() const {
312 return m_bounds;
313 }
314
315 const std::vector<double> & Timeseries::times() const {
316 return m_time;
317 }
318
319 const std::vector<double> & Timeseries::time_bounds() const {
320 return m_time_bounds;
321 }
322
323 const std::vector<double> & Timeseries::values() const {
324 return m_values;
325 }
326
327 } // end of namespace pism