PISM, A Parallel Ice Sheet Model
ScalarForcing.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2019, 2020, 2021 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 "ScalarForcing.hh"
21 
22 #include <algorithm> // std::min_element(), std::max_element()
23 #include <cassert> // assert()
24 #include <cmath> // std::floor()
25 
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/util/Context.hh"
28 #include "pism/util/Time.hh"
29 #include "pism/util/error_handling.hh"
30 #include "pism/util/Logger.hh"
31 #include "pism/util/io/File.hh"
32 #include "pism/util/io/io_helpers.hh"
33 #include "pism/util/VariableMetadata.hh"
34 
35 namespace pism {
36 
37 //! \brief Report the range of a time-series stored in `data`.
38 static void report_range(const std::vector<double> &data,
39  const units::System::Ptr &unit_system,
40  const VariableMetadata &metadata,
41  const Logger &log) {
42  // min_element and max_element return iterators; "*" is used to get
43  // the value corresponding to this iterator
44  double min = *std::min_element(data.begin(), data.end());
45  double max = *std::max_element(data.begin(), data.end());
46 
47  units::Converter c(unit_system,
48  metadata.get_string("units"),
49  metadata.get_string("glaciological_units"));
50  min = c(min);
51  max = c(max);
52 
53  std::string spacer(metadata.get_name().size(), ' ');
54 
55  log.message(2,
56  " FOUND %s / %-60s\n"
57  " %s \\ min,max = %9.3f,%9.3f (%s)\n",
58  metadata.get_name().c_str(),
59  metadata.get_string("long_name").c_str(), spacer.c_str(), min, max,
60  metadata.get_string("glaciological_units").c_str());
61 }
62 
64  const std::string &filename,
65  const std::string &variable_name,
66  const std::string &units,
67  const std::string &glaciological_units,
68  const std::string &long_name,
69  bool periodic) {
70  try {
71  auto unit_system = ctx.unit_system();
72 
73  VariableMetadata variable(variable_name, unit_system);
74 
75  variable["units"] = units;
76  variable["glaciological_units"] = glaciological_units;
77  variable["long_name"] = long_name;
78 
79  double forcing_t0{};
80  double forcing_t1{};
81 
82  // Read data from a NetCDF file
83  {
84  ctx.log()->message(2,
85  " reading %s (%s) from file '%s'...\n",
86  long_name.c_str(), variable_name.c_str(), filename.c_str());
87 
88  File file(ctx.com(), filename, PISM_NETCDF3, PISM_READONLY);
89 
90  // Read forcing data. The read_timeseries() call will ensure that variable_name is a
91  // scalar variable.
92  std::vector<double> data{};
93  io::read_timeseries(file, variable, *ctx.log(), data);
94 
95  // The following line relies on the fact that this is a scalar variable.
96  std::string time_name = file.dimensions(variable_name)[0];
97 
98  std::vector<double> times{};
99  std::vector<double> bounds{};
100  io::read_time_info(*ctx.log(), unit_system,
101  file, time_name, ctx.time()->units_string(),
102  times, bounds);
103  size_t N = times.size();
104 
105  // Compute values used to extend data read from file
106  //
107  // Initialize using values corresponding to constant extrapolation:
108  double v0 = data.front();
109  double v1 = data.back();
110  if (periodic) {
111  double b0 = bounds.front();
112  double b1 = bounds.back();
113 
114  // compute the value at the beginning and end of the period:
115  {
116  double dt1 = times.front() - b0;
117  double dt2 = b1 - times.back();
118 
119  // interpolation factor
120  double alpha = dt1 + dt2 > 0 ? dt2 / (dt1 + dt2) : 0.0;
121 
122  v0 = (1.0 - alpha) * data.back() + alpha * data.front();
123  v1 = v0;
124  }
125 
126  {
127  m_period = b1 - b0;
128  m_period_start = b0;
129 
130  assert(m_period > 0.0);
131  }
132  }
133 
134  // Note: this should take care of file with one record as well.
135  m_times.clear();
136  m_values.clear();
137  {
138  if (bounds.front() < times.front()) {
139  m_times.emplace_back(bounds.front());
140  m_values.emplace_back(v0);
141  }
142  for (size_t k = 0; k < N; ++k) {
143  m_times.emplace_back(times[k]);
144  m_values.emplace_back(data[k]);
145  }
146  if (bounds.back() > times.back()) {
147  m_times.emplace_back(bounds.back());
148  m_values.emplace_back(v1);
149  }
150  }
151 
152  // Save the time interval covered by this forcing.
153  forcing_t0 = bounds.front();
154  forcing_t1 = bounds.back();
155  } // end of the block initializing data
156 
157  // validate resulting times
158  if (not is_increasing(m_times)) {
160  "m_times have to be strictly increasing (this is a bug)");
161  }
162 
163  report_range(m_values, unit_system, variable, *ctx.log());
164 
165  // Set up interpolation.
166  {
167  m_acc = gsl_interp_accel_alloc();
168  m_spline = gsl_spline_alloc(gsl_interp_linear, m_times.size());
169  gsl_spline_init(m_spline, m_times.data(), m_values.data(), m_times.size());
170  }
171 
172  bool extrapolate = ctx.config()->get_flag("input.forcing.time_extrapolation");
173 
174  if (not (extrapolate or periodic)) {
175  check_forcing_duration(*ctx.time(), forcing_t0, forcing_t1);
176  }
177 
178  } catch (RuntimeError &e) {
179  e.add_context("reading %s (%s) from '%s'",
180  long_name.c_str(), variable_name.c_str(), filename.c_str());
181  throw;
182  }
183 }
184 
186  const std::string &prefix,
187  const std::string &variable_name,
188  const std::string &units,
189  const std::string &glaciological_units,
190  const std::string &long_name)
191  : m_period(0.0),
192  m_period_start(0.0),
193  m_acc(nullptr),
194  m_spline(nullptr) {
195 
196  Config::ConstPtr config = ctx.config();
197 
198  auto filename = config->get_string(prefix + ".file");
199  bool periodic = config->get_flag(prefix + ".periodic");
200 
201  if (filename.empty()) {
203  "%s.file is required", prefix.c_str());
204  }
205 
206  initialize(ctx, filename, variable_name, units, glaciological_units,
207  long_name, periodic);
208 }
209 
211  const std::string &filename,
212  const std::string &variable_name,
213  const std::string &units,
214  const std::string &glaciological_units,
215  const std::string &long_name,
216  bool periodic)
217  : m_period(0.0),
218  m_period_start(0.0),
219  m_acc(nullptr),
220  m_spline(nullptr) {
221 
222  initialize(ctx, filename, variable_name, units, glaciological_units,
223  long_name, periodic);
224 }
225 
227  if (m_spline != nullptr) {
228  gsl_spline_free(m_spline);
229  }
230  if (m_acc != nullptr) {
231  gsl_interp_accel_free(m_acc);
232  }
233 }
234 
235 double ScalarForcing::value(double t) const {
236  double T = t;
237 
238  if (m_period > 0.0) {
239  // compute time since the period start
240  T -= m_period_start;
241 
242  // remove an integer number of periods
243  T -= std::floor(T / m_period) * m_period;
244 
245  // add the period start back
246  T += m_period_start;
247  }
248 
249  if (T > m_times.back()) {
250  return m_values.back();
251  }
252 
253  if (T < m_times.front()) {
254  return m_values.front();
255  }
256 
257  return gsl_spline_eval(m_spline, T, m_acc);
258 }
259 
260 // Integrate from a to b, interpreting data as *not* periodic
261 double ScalarForcing::integral(double a, double b) const {
262  assert(b >= a);
263 
264  double dt = b - a;
265 
266  double t0 = m_times.front();
267  double t1 = m_times.back();
268  double v0 = m_values[0];
269  double v1 = m_values.back();
270 
271  // both points are to the left of [t0, t1]
272  if (b <= t0) {
273  return v0 * (b - a);
274  }
275 
276  // both points are to the right of [t0, t1]
277  if (a >= t1) {
278  return v1 * (b - a);
279  }
280 
281  // a is to the left of [t0, t1]
282  if (a < t0) {
283  return v0 * (t0 - a) + integral(t0, b);
284  }
285 
286  // b is to the right of [t0, t1]
287  if (b > t1) {
288  return integral(a, t1) + v1 * (b - t1);
289  }
290 
291  // both points are inside [t0, t1]
292  size_t ai = gsl_interp_bsearch(m_times.data(), a, 0, m_times.size() - 1);
293  size_t bi = gsl_interp_bsearch(m_times.data(), b, 0, m_times.size() - 1);
294 
295  // gsl_interp_bsearch() above returns the index i of the array ‘m_times’ such that
296  // ‘m_times[i] <= t < m_times[i+1]’. The index is searched for in the
297  // range [0, m_times.size() - 1] (inclusive).
298 
299  double v_a = gsl_spline_eval(m_spline, a, m_acc);
300  double v_b = gsl_spline_eval(m_spline, b, m_acc);
301 
302  if (ai == bi) {
303  return 0.5 * (v_a + v_b) * dt;
304  }
305 
306  double result = 0.0;
307  // integrate from a to the data point just to its right
308  result += 0.5 * (v_a + m_values[ai + 1]) * (m_times[ai + 1] - a);
309 
310  // integrate over (possibly zero) intervals between a and b
311  for (size_t k = ai + 1; k < bi; ++k) {
312  result += 0.5 * (m_values[k] + m_values[k + 1]) * (m_times[k + 1] - m_times[k]);
313  }
314 
315  result += 0.5 * (m_values[bi] + v_b) * (b - m_times[bi]);
316 
317  return result;
318 }
319 
320 double ScalarForcing::average(double t, double dt) const {
321  if (dt == 0.0) {
322  return value(t);
323  }
324 
325  if (dt < 0.0) {
326  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "negative interval length");
327  }
328 
329  if (m_period <= 0.0) {
330  // regular case
331  return integral(t, t + dt) / dt;
332  }
333 
334  // periodic case
335  {
336  double a = t;
337  double b = t + dt;
338  double t0 = m_period_start;
339  double P = m_period;
340 
341  double N = std::floor((a - t0) / P);
342  double M = std::floor((b - t0) / P);
343  double delta = a - t0 - P * N;
344  double gamma = b - t0 - P * M;
345 
346  double N_periods = M - (N + 1);
347 
348  if (N_periods >= 0.0) {
349  return (integral(t0 + delta, t0 + P) +
350  N_periods * integral(t0, t0 + P) +
351  integral(t0, t0 + gamma)) / dt;
352  }
353 
354  return integral(t0 + delta, t0 + gamma) / dt;
355  }
356 }
357 
358 
359 } // end of namespace pism
std::shared_ptr< const Config > ConstPtr
MPI_Comm com() const
Definition: Context.cc:82
ConstConfigPtr config() const
Definition: Context.cc:106
UnitsSystemPtr unit_system() const
Definition: Context.cc:98
ConstLoggerPtr log() const
Definition: Context.cc:130
ConstTimePtr time() const
Definition: Context.cc:118
std::vector< std::string > dimensions(const std::string &variable_name) const
Definition: File.cc:428
High-level PISM I/O class.
Definition: File.hh:51
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
Definition: Logger.cc:50
A basic logging class.
Definition: Logger.hh:40
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
gsl_spline * m_spline
double average(double t, double dt) const
double value(double t) const
std::vector< double > m_times
std::vector< double > m_values
ScalarForcing(const Context &ctx, const std::string &option_prefix, const std::string &variable_name, const std::string &units, const std::string &glaciological_units, const std::string &long_name)
void initialize(const Context &ctx, const std::string &filename, const std::string &variable_name, const std::string &units, const std::string &glaciological_units, const std::string &long_name, bool periodic)
gsl_interp_accel * m_acc
double integral(double a, double b) const
std::string get_string(const std::string &name) const
Get a string attribute.
std::string get_name() const
std::shared_ptr< System > Ptr
Definition: Units.hh:47
#define PISM_ERROR_LOCATION
static const double b0
Definition: exactTestL.cc:41
void read_timeseries(const File &file, const VariableMetadata &metadata, const Logger &log, std::vector< double > &data)
Read a time-series variable from a NetCDF file to a vector of doubles.
Definition: io_helpers.cc:1049
void read_time_info(const Logger &log, std::shared_ptr< units::System > unit_system, const File &file, const std::string &time_name, const std::string &time_units, std::vector< double > &times, std::vector< double > &bounds)
Definition: io_helpers.cc:1286
bool is_increasing(const std::vector< double > &a)
Checks if a vector of doubles is strictly increasing.
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double v0
Definition: exactTestP.cc:53
static const double k
Definition: exactTestP.cc:45
static void report_range(const std::vector< double > &data, const units::System::Ptr &unit_system, const VariableMetadata &metadata, const Logger &log)
Report the range of a time-series stored in data.
void check_forcing_duration(const Time &time, double forcing_start, double forcing_end)
Definition: Time.cc:932
@ PISM_NETCDF3
Definition: IO_Flags.hh:41
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49