ticeModelVec2T.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
---
ticeModelVec2T.cc (18844B)
---
1 // Copyright (C) 2009--2019 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 <petsc.h>
20 #include <cassert>
21
22 #include "iceModelVec2T.hh"
23 #include "pism/util/io/File.hh"
24 #include "pism_utilities.hh"
25 #include "Time.hh"
26 #include "IceGrid.hh"
27 #include "ConfigInterface.hh"
28
29 #include "error_handling.hh"
30 #include "io/io_helpers.hh"
31 #include "pism/util/Logger.hh"
32 #include "pism/util/interpolation.hh"
33
34 namespace pism {
35
36
37 /*!
38 * Allocate an instance that will be used to load and use a forcing field from a file.
39 *
40 * Checks the number of records in a file and allocates storage accordingly.
41 *
42 * If `periodic` is true, allocate enough storage to hold all the records, otherwise
43 * allocate storage for at most `max_buffer_size` records.
44 *
45 * @param[in] grid computational grid
46 * @param[in] file input file
47 * @param[in] short_name variable name in `file`
48 * @param[in] standard_name standard name (if available); leave blank to ignore
49 * @param[in] max_buffer_size maximum buffer size for non-periodic fields
50 * @param[in] evaluations_per_year number of evaluations per year to use when averaging
51 * @param[in] periodic true if this forcing field should be interpreted as periodic
52 */
53 IceModelVec2T::Ptr IceModelVec2T::ForcingField(IceGrid::ConstPtr grid,
54 const File &file,
55 const std::string &short_name,
56 const std::string &standard_name,
57 int max_buffer_size,
58 int evaluations_per_year,
59 bool periodic,
60 InterpolationType interpolation_type) {
61
62 int n_records = file.nrecords(short_name, standard_name,
63 grid->ctx()->unit_system());
64
65 if (not periodic) {
66 n_records = std::min(n_records, max_buffer_size);
67 }
68 // In the periodic case we try to keep all the records in RAM.
69
70 // Allocate storage for one record if the variable was not found. This is needed to be
71 // able to cheaply allocate and then discard an "-atmosphere given" model
72 // (atmosphere::Given) when "-surface given" (Given) is selected.
73 n_records = std::max(n_records, 1);
74
75 // LCOV_EXCL_START
76 if (n_records > IceGrid::max_dm_dof) {
77 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
78 "cannot allocate storage for %d records of %s (%s)"
79 " (exceeds the maximum of %d)",
80 n_records, short_name.c_str(), standard_name.c_str(),
81 IceGrid::max_dm_dof);
82 }
83 // LCOV_EXCL_STOP
84
85 if (periodic and interpolation_type == LINEAR) {
86 interpolation_type = LINEAR_PERIODIC;
87 }
88
89 return IceModelVec2T::Ptr(new IceModelVec2T(grid, short_name, n_records,
90 evaluations_per_year, interpolation_type));
91 }
92
93
94 IceModelVec2T::IceModelVec2T(IceGrid::ConstPtr grid, const std::string &short_name,
95 unsigned int n_records,
96 unsigned int n_evaluations_per_year,
97 InterpolationType interpolation_type)
98 : IceModelVec2S(grid, short_name, WITHOUT_GHOSTS, 1),
99 m_array3(nullptr),
100 m_n_records(n_records),
101 m_N(0),
102 m_n_evaluations_per_year(n_evaluations_per_year),
103 m_first(-1),
104 m_interp_type(interpolation_type),
105 m_period(0),
106 m_reference_time(0.0)
107 {
108 m_report_range = false;
109
110 if (not (m_interp_type == PIECEWISE_CONSTANT or
111 m_interp_type == LINEAR or
112 m_interp_type == LINEAR_PERIODIC)) {
113 throw RuntimeError(PISM_ERROR_LOCATION, "unsupported interpolation type");
114 }
115
116 // LCOV_EXCL_START
117 if (n_records > IceGrid::max_dm_dof) {
118 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
119 "cannot allocate storage for %d records of %s"
120 " (exceeds the maximum of %d)",
121 n_records, short_name.c_str(), IceGrid::max_dm_dof);
122 }
123 // LCOV_EXCL_STOP
124
125 // initialize the m_da3 member:
126 m_da3 = m_grid->get_dm(n_records, this->m_da_stencil_width);
127
128 // allocate the 3D Vec:
129 PetscErrorCode ierr = DMCreateGlobalVector(*m_da3, m_v3.rawptr());
130 PISM_CHK(ierr, "DMCreateGlobalVector");
131 }
132
133 IceModelVec2T::~IceModelVec2T() {
134 // empty
135 }
136
137 unsigned int IceModelVec2T::n_records() {
138 return m_n_records;
139 }
140
141 double*** IceModelVec2T::get_array3() {
142 begin_access();
143 return reinterpret_cast<double***>(m_array3);
144 }
145
146 void IceModelVec2T::begin_access() const {
147 if (m_access_counter == 0) {
148 PetscErrorCode ierr = DMDAVecGetArrayDOF(*m_da3, m_v3, &m_array3);
149 PISM_CHK(ierr, "DMDAVecGetArrayDOF");
150 }
151
152 // this call will increment the m_access_counter
153 IceModelVec2S::begin_access();
154
155 }
156
157 void IceModelVec2T::end_access() const {
158 // this call will decrement the m_access_counter
159 IceModelVec2S::end_access();
160
161 if (m_access_counter == 0) {
162 PetscErrorCode ierr = DMDAVecRestoreArrayDOF(*m_da3, m_v3, &m_array3);
163 PISM_CHK(ierr, "DMDAVecRestoreArrayDOF");
164 m_array3 = NULL;
165 }
166 }
167
168 void IceModelVec2T::init(const std::string &fname, unsigned int period, double reference_time) {
169
170 const Logger &log = *m_grid->ctx()->log();
171
172 m_filename = fname;
173 m_period = period;
174 m_reference_time = reference_time;
175
176 // We find the variable in the input file and
177 // try to find the corresponding time dimension.
178
179 File file(m_grid->com, m_filename, PISM_GUESS, PISM_READONLY);
180 auto var = file.find_variable(m_metadata[0].get_name(), m_metadata[0].get_string("standard_name"));
181 if (not var.exists) {
182 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "can't find %s (%s) in %s.",
183 m_metadata[0].get_string("long_name").c_str(),
184 m_metadata[0].get_name().c_str(),
185 m_filename.c_str());
186 }
187
188 auto time_name = io::time_dimension(m_grid->ctx()->unit_system(),
189 file, var.name);
190
191 if (not time_name.empty()) {
192 // we're found the time dimension
193 TimeseriesMetadata time_dimension(time_name, time_name, m_grid->ctx()->unit_system());
194
195 auto time_units = m_grid->ctx()->time()->units_string();
196 time_dimension.set_string("units", time_units);
197
198 io::read_timeseries(file, time_dimension,
199 *m_grid->ctx()->time(), log, m_time);
200
201 std::string bounds_name = file.read_text_attribute(time_name, "bounds");
202
203 if (m_time.size() > 1) {
204
205 if (m_interp_type == PIECEWISE_CONSTANT) {
206 if (bounds_name.empty()) {
207 // no time bounds attribute
208 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
209 "Variable '%s' does not have the time_bounds attribute.\n"
210 "Cannot use time-dependent forcing data '%s' (%s) without time bounds.",
211 time_name.c_str(), m_metadata[0].get_string("long_name").c_str(),
212 m_metadata[0].get_name().c_str());
213 }
214
215 // read time bounds data from a file
216 TimeBoundsMetadata tb(bounds_name, time_name, m_grid->ctx()->unit_system());
217 tb.set_string("units", time_units);
218
219 io::read_time_bounds(file, tb, *m_grid->ctx()->time(),
220 log, m_time_bounds);
221
222 // time bounds data overrides the time variable: we make t[j] be the
223 // left end-point of the j-th interval
224 for (unsigned int k = 0; k < m_time.size(); ++k) {
225 m_time[k] = m_time_bounds[2*k + 0];
226 }
227 } else {
228 // fake time step length used to generate the right end point of the last interval
229 // TODO: figure out if there is a better way to do this.
230 double dt = 1.0;
231 size_t N = m_time.size();
232 m_time_bounds.resize(2 * N);
233 for (size_t k = 0; k < N; ++k) {
234 m_time_bounds[2 * k + 0] = m_time[k];
235 m_time_bounds[2 * k + 1] = k + 1 < N ? m_time[k + 1] : m_time[k] + dt;
236 }
237 }
238
239 } else {
240 // only one time record; set fake time bounds:
241 m_time_bounds = {m_time[0] - 1.0, m_time[0] + 1};
242 }
243
244 } else {
245 // no time dimension; assume that we have only one record and set the time
246 // to 0
247 m_time = {0.0};
248
249 // set fake time bounds:
250 m_time_bounds = {-1.0, 1.0};
251 }
252
253 if (not is_increasing(m_time)) {
254 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
255 "times have to be strictly increasing (read from '%s').",
256 m_filename.c_str());
257 }
258
259 if (m_period != 0) {
260 if ((size_t)m_n_records < m_time.size()) {
261 throw RuntimeError(PISM_ERROR_LOCATION,
262 "buffer has to be big enough to hold all records of periodic data");
263 }
264
265 // read periodic data right away (we need to hold it all in memory anyway)
266 update(0);
267 }
268 }
269
270 //! Initialize as constant in time and space
271 void IceModelVec2T::init_constant(double value) {
272
273 // set constant value everywhere
274 set(value);
275 set_record(0);
276
277 // set the time to zero
278 m_time = {0.0};
279 m_N = 1;
280 m_first = 0;
281
282 // set fake time bounds:
283 m_time_bounds = {-1.0, 1.0};
284 }
285
286 //! Read some data to make sure that the interval (t, t + dt) is covered.
287 void IceModelVec2T::update(double t, double dt) {
288
289 if (m_filename.empty()) {
290 // We are not reading data from a file.
291 return;
292 }
293
294 if (m_time_bounds.size() == 0) {
295 update(0);
296 return;
297 }
298
299 if (m_period != 0) {
300 // we read all data in IceModelVec2T::init() (see above)
301 return;
302 }
303
304 if (m_N > 0) {
305 unsigned int last = m_first + (m_N - 1);
306
307 // find the interval covered by data held in memory:
308 double
309 t0 = m_time_bounds[m_first * 2],
310 t1 = m_time_bounds[last * 2 + 1];
311
312 // just return if we have all the data we need:
313 if (t >= t0 and t + dt <= t1) {
314 return;
315 }
316 }
317
318 Interpolation I(m_interp_type, m_time, {t, t + dt});
319
320 unsigned int
321 first = I.left(0),
322 last = I.right(1),
323 N = last - first + 1;
324
325 // check if all the records necessary to cover this interval fit in the
326 // buffer:
327 if (N > m_n_records) {
328 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
329 "cannot read %d records of %s (buffer size: %d)",
330 N, m_name.c_str(), m_n_records);
331 }
332
333 update(first);
334 }
335
336 //! Update by reading at most n_records records from the file.
337 void IceModelVec2T::update(unsigned int start) {
338
339 unsigned int time_size = (int)m_time.size();
340
341 if (start >= time_size) {
342 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
343 "IceModelVec2T::update(int start): start = %d is invalid", start);
344 }
345
346 unsigned int missing = std::min(m_n_records, time_size - start);
347
348 if (start == static_cast<unsigned int>(m_first)) {
349 // nothing to do
350 return;
351 }
352
353 int kept = 0;
354 if (m_first >= 0) {
355 unsigned int last = m_first + (m_N - 1);
356 if ((m_N > 0) && (start >= (unsigned int)m_first) && (start <= last)) {
357 int discarded = start - m_first;
358 kept = last - start + 1;
359 discard(discarded);
360 missing -= kept;
361 start += kept;
362 m_first += discarded;
363 } else {
364 m_first = start;
365 }
366 } else {
367 m_first = start;
368 }
369
370 if (missing <= 0) {
371 return;
372 }
373
374 m_N = kept + missing;
375
376 Time::ConstPtr t = m_grid->ctx()->time();
377
378 Logger::ConstPtr log = m_grid->ctx()->log();
379 if (this->n_records() > 1) {
380 log->message(4,
381 " reading \"%s\" into buffer\n"
382 " (short_name = %s): %d records, time intervals (%s, %s) through (%s, %s)...\n",
383 metadata().get_string("long_name").c_str(), m_name.c_str(), missing,
384 t->date(m_time_bounds[start*2]).c_str(),
385 t->date(m_time_bounds[start*2 + 1]).c_str(),
386 t->date(m_time_bounds[(start + missing - 1)*2]).c_str(),
387 t->date(m_time_bounds[(start + missing - 1)*2 + 1]).c_str());
388 m_report_range = false;
389 } else {
390 m_report_range = true;
391 }
392
393 File file(m_grid->com, m_filename, PISM_GUESS, PISM_READONLY);
394
395 const bool allow_extrapolation = m_grid->ctx()->config()->get_flag("grid.allow_extrapolation");
396
397 for (unsigned int j = 0; j < missing; ++j) {
398 {
399 petsc::VecArray tmp_array(m_v);
400 io::regrid_spatial_variable(m_metadata[0], *m_grid, file, start + j, CRITICAL,
401 m_report_range, allow_extrapolation,
402 0.0, m_interpolation_type, tmp_array.get());
403 }
404
405 m_grid->ctx()->log()->message(5, " %s: reading entry #%02d, year %s...\n",
406 m_name.c_str(),
407 start + j,
408 t->date(m_time[start + j]).c_str());
409
410 set_record(kept + j);
411 }
412 }
413
414 //! Discard the first N records, shifting the rest of them towards the "beginning".
415 void IceModelVec2T::discard(int number) {
416
417 if (number == 0) {
418 return;
419 }
420
421 m_N -= number;
422
423 double ***a3 = get_array3();
424 for (Points p(*m_grid); p; p.next()) {
425 const int i = p.i(), j = p.j();
426
427 for (unsigned int k = 0; k < m_N; ++k) {
428 a3[j][i][k] = a3[j][i][k + number];
429 }
430 }
431 end_access();
432 }
433
434 //! Sets the record number n to the contents of the (internal) Vec v.
435 void IceModelVec2T::set_record(int n) {
436
437 double **a2 = get_array();
438 double ***a3 = get_array3();
439 for (Points p(*m_grid); p; p.next()) {
440 const int i = p.i(), j = p.j();
441 a3[j][i][n] = a2[j][i];
442 }
443 end_access();
444 end_access();
445 }
446
447 //! Sets the (internal) Vec v to the contents of the nth record.
448 void IceModelVec2T::get_record(int n) {
449
450 double **a2 = get_array();
451 double ***a3 = get_array3();
452 for (Points p(*m_grid); p; p.next()) {
453 const int i = p.i(), j = p.j();
454 a2[j][i] = a3[j][i][n];
455 }
456 end_access();
457 end_access();
458 }
459
460 //! @brief Given the time t determines the maximum possible time-step this IceModelVec2T
461 //! allows.
462 MaxTimestep IceModelVec2T::max_timestep(double t) const {
463 // only allow going to the next record
464
465 // find the index k such that m_time[k] <= x < m_time[k + 1]
466 size_t k = gsl_interp_bsearch(m_time.data(), t, 0, m_time.size());
467
468 // end of the corresponding interval
469 double
470 t_next = m_time_bounds[2 * k + 1],
471 dt = std::max(t_next - t, 0.0);
472
473 if (dt > 1.0) { // never take time-steps shorter than 1 second
474 return MaxTimestep(dt);
475 } else if (k + 1 < m_time.size()) {
476 dt = m_time_bounds[2 * (k + 1) + 1] - m_time_bounds[2 * (k + 1)];
477 return MaxTimestep(dt);
478 } else {
479 return MaxTimestep();
480 }
481 }
482
483 /*
484 * \brief Use piecewise-constant interpolation to initialize
485 * IceModelVec2T with the value at time `t`.
486 *
487 * \note This method does not check if an update() call is necessary!
488 *
489 * @param[in] t requested time
490 *
491 */
492 void IceModelVec2T::interp(double t) {
493
494 init_interpolation({t});
495
496 get_record(m_interp->left(0));
497 }
498
499
500 /**
501 * Compute the average value over the time interval `[t, t + dt]`.
502 *
503 * @param t start of the time interval, in seconds
504 * @param dt length of the time interval, in seconds
505 *
506 */
507 void IceModelVec2T::average(double t, double dt) {
508
509 double dt_years = units::convert(m_grid->ctx()->unit_system(),
510 dt, "seconds", "years"); // *not* time->year(dt)
511
512 // if only one record, nothing to do
513 if (m_time.size() == 1) {
514 return;
515 }
516
517 // Determine the number of small time-steps to use for averaging:
518 int M = (int) ceil(m_n_evaluations_per_year * (dt_years));
519 if (M < 1) {
520 M = 1;
521 }
522
523 std::vector<double> ts(M);
524 double ts_dt = dt / M;
525 for (int k = 0; k < M; k++) {
526 ts[k] = t + k * ts_dt;
527 }
528
529 init_interpolation(ts);
530
531 double **a2 = get_array(); // calls begin_access()
532 for (Points p(*m_grid); p; p.next()) {
533 const int i = p.i(), j = p.j();
534 a2[j][i] = average(i, j);
535 }
536 end_access();
537 }
538
539 /**
540 * \brief Compute weights for the piecewise-constant interpolation.
541 * This is used *both* for time-series and "snapshots".
542 *
543 * @param ts requested times, in seconds
544 *
545 */
546 void IceModelVec2T::init_interpolation(const std::vector<double> &ts) {
547
548 assert(m_first >= 0);
549
550 auto time = m_grid->ctx()->time();
551
552 // Compute "periodized" times if necessary.
553 std::vector<double> times_requested(ts.size());
554 if (m_period != 0) {
555 for (unsigned int k = 0; k < ts.size(); ++k) {
556 times_requested[k] = time->mod(ts[k] - m_reference_time, m_period);
557 }
558 } else {
559 times_requested = ts;
560 }
561
562 m_interp.reset(new Interpolation(m_interp_type, &m_time[m_first], m_N,
563 times_requested.data(), times_requested.size(),
564 time->years_to_seconds(m_period)));
565 }
566
567 /**
568 * \brief Compute values of the time-series using precomputed indices
569 * (and piecewise-constant interpolation).
570 *
571 * @param i,j map-plane grid point
572 * @param result pointer to an allocated array of `weights.size()` `double`
573 *
574 */
575 void IceModelVec2T::interp(int i, int j, std::vector<double> &result) {
576 double ***a3 = (double***) m_array3;
577
578 result.resize(m_interp->alpha().size());
579
580 m_interp->interpolate(a3[j][i], result.data());
581 }
582
583 //! \brief Finds the average value at i,j over the interval (t, t +
584 //! dt) using the rectangle rule.
585 /*!
586 Can (and should) be optimized. Later, though.
587 */
588 double IceModelVec2T::average(int i, int j) {
589 unsigned int M = m_interp->alpha().size();
590 double result = 0.0;
591
592 if (m_N == 1) {
593 double ***a3 = (double***) m_array3;
594 result = a3[j][i][0];
595 } else {
596 std::vector<double> values(M);
597
598 interp(i, j, values);
599
600 // rectangular rule (uses the fact that points are equally-spaced
601 // in time)
602 result = 0;
603 for (unsigned int k = 0; k < M; ++k) {
604 result += values[k];
605 }
606 result /= (double)M;
607 }
608 return result;
609 }
610
611
612
613 } // end of namespace pism