PISM, A Parallel Ice Sheet Model  stable v2.1.1 committed by Constantine Khrulev on 2024-12-04 13:36:58 -0900
Diagnostic.cc
Go to the documentation of this file.
1 /* Copyright (C) 2015, 2016, 2017, 2019, 2020, 2021, 2022, 2023 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 #include <cmath>
20 
21 #include "pism/util/Diagnostic.hh"
22 #include "pism/util/Time.hh"
23 #include "pism/util/error_handling.hh"
24 #include "pism/util/io/io_helpers.hh"
25 #include "pism/util/Logger.hh"
26 #include "pism/util/pism_utilities.hh"
27 #include "pism/util/Context.hh"
28 
29 namespace pism {
30 
31 Diagnostic::Diagnostic(std::shared_ptr<const Grid> grid)
32  : m_grid(grid),
33  m_sys(grid->ctx()->unit_system()),
34  m_config(grid->ctx()->config()),
35  m_fill_value(m_config->get_number("output.fill_value")) {
36  // empty
37 }
38 
39 void Diagnostic::update(double dt) {
40  this->update_impl(dt);
41 }
42 
43 void Diagnostic::update_impl(double dt) {
44  (void) dt;
45  // empty
46 }
47 
49  this->reset_impl();
50 }
51 
53  // empty
54 }
55 
56 /*!
57  * Convert from external (output) units to internal units.
58  */
59 double Diagnostic::to_internal(double x) const {
60  std::string
61  out = m_vars.at(0)["output_units"],
62  in = m_vars.at(0)["units"];
63  return convert(m_sys, x, out, in);
64 }
65 
66 /*!
67  * Convert from internal to external (output) units.
68  */
69 double Diagnostic::to_external(double x) const {
70  std::string
71  out = m_vars.at(0)["output_units"],
72  in = m_vars.at(0)["units"];
73  return convert(m_sys, x, in, out);
74 }
75 
76 //! Get the number of NetCDF variables corresponding to a diagnostic quantity.
77 unsigned int Diagnostic::n_variables() const {
78  return m_vars.size();
79 }
80 
81 void Diagnostic::init(const File &input, unsigned int time) {
82  this->init_impl(input, time);
83 }
84 
85 void Diagnostic::define_state(const File &output) const {
86  this->define_state_impl(output);
87 }
88 
89 void Diagnostic::write_state(const File &output) const {
90  this->write_state_impl(output);
91 }
92 
93 void Diagnostic::init_impl(const File &input, unsigned int time) {
94  (void) input;
95  (void) time;
96  // empty
97 }
98 
99 void Diagnostic::define_state_impl(const File &output) const {
100  (void) output;
101  // empty
102 }
103 
104 void Diagnostic::write_state_impl(const File &output) const {
105  (void) output;
106  // empty
107 }
108 
109 //! Get a metadata object corresponding to variable number N.
111  if (N >= m_vars.size()) {
113  "variable metadata index %d is out of bounds",
114  N);
115  }
116 
117  return m_vars[N];
118 }
119 
120 void Diagnostic::define(const File &file, io::Type default_type) const {
121  this->define_impl(file, default_type);
122 }
123 
124 //! Define NetCDF variables corresponding to a diagnostic quantity.
125 void Diagnostic::define_impl(const File &file, io::Type default_type) const {
126  for (const auto &v : m_vars) {
127  io::define_spatial_variable(v, *m_grid, file, default_type);
128  }
129 }
130 
131 std::shared_ptr<array::Array> Diagnostic::compute() const {
132  std::vector<std::string> names;
133  for (const auto &v : m_vars) {
134  names.push_back(v.get_name());
135  }
136  std::string all_names = join(names, ",");
137 
138  m_grid->ctx()->log()->message(3, "- Computing %s...\n", all_names.c_str());
139  auto result = this->compute_impl();
140  m_grid->ctx()->log()->message(3, "- Done computing %s.\n", all_names.c_str());
141 
142  return result;
143 }
144 
145 TSDiagnostic::TSDiagnostic(std::shared_ptr<const Grid> grid, const std::string &name)
146  : m_grid(grid),
147  m_config(grid->ctx()->config()),
148  m_sys(grid->ctx()->unit_system()),
149  m_time_name(grid->ctx()->config()->get_string("time.dimension_name")),
150  m_variable(name, m_sys),
151  m_dimension(m_time_name, m_sys),
152  m_time_bounds(m_time_name + "_bounds", m_sys) {
153 
154  m_current_time = 0;
155  m_start = 0;
156 
157  m_buffer_size = static_cast<size_t>(m_config->get_number("output.timeseries.buffer_size"));
158 
159  m_variable["ancillary_variables"] = name + "_aux";
160 
161  m_dimension["calendar"] = m_grid->ctx()->time()->calendar();
162  m_dimension["long_name"] = "time";
163  m_dimension["axis"] = "T";
164  m_dimension["units"] = m_grid->ctx()->time()->units_string();
165 }
166 
168  flush();
169 }
170 
171 void TSDiagnostic::set_units(const std::string &units,
172  const std::string &output_units) {
173  m_variable["units"] = units;
174 
175  if (not m_config->get_flag("output.use_MKS")) {
176  m_variable["output_units"] = output_units;
177  }
178 }
179 
180 TSSnapshotDiagnostic::TSSnapshotDiagnostic(std::shared_ptr<const Grid> grid, const std::string &name)
181  : TSDiagnostic(grid, name) {
182  // empty
183 }
184 
185 TSRateDiagnostic::TSRateDiagnostic(std::shared_ptr<const Grid> grid, const std::string &name)
186  : TSDiagnostic(grid, name), m_accumulator(0.0), m_v_previous(0.0), m_v_previous_set(false) {
187  // empty
188 }
189 
190 TSFluxDiagnostic::TSFluxDiagnostic(std::shared_ptr<const Grid> grid, const std::string &name)
191  : TSRateDiagnostic(grid, name) {
192  // empty
193 }
194 
195 void TSSnapshotDiagnostic::evaluate(double t0, double t1, double v) {
196 
197  // skip times before the beginning of this time step
198  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] < t0) {
199  m_current_time += 1;
200  }
201 
202  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] <= t1) {
203  const unsigned int k = m_current_time;
204  m_current_time += 1;
205 
206  // skip the first time: it defines the beginning of a reporting interval
207  if (k == 0) {
208  continue;
209  }
210 
211  const double t_s = (*m_requested_times)[k - 1];
212  const double t_e = (*m_requested_times)[k];
213 
214  // store computed data in the buffer
215  {
216  m_time.push_back(t_e);
217  m_values.push_back(v);
218  m_bounds.push_back(t_s);
219  m_bounds.push_back(t_e);
220  }
221  }
222 }
223 
224 void TSRateDiagnostic::evaluate(double t0, double t1, double change) {
225  static const double epsilon = 1e-4; // seconds
226  assert(t1 > t0);
227 
228  // skip times before and including the beginning of this time step
229  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] <= t0) {
230  m_current_time += 1;
231  }
232 
233  // number of requested times in this time step
234  unsigned int N = 0;
235 
236  // loop through requested times that are within this time step
237  while (m_current_time < m_requested_times->size() and (*m_requested_times)[m_current_time] <= t1) {
238  const unsigned int k = m_current_time;
239  m_current_time += 1;
240 
241  N += 1;
242 
243  // skip the first time: it defines the beginning of a reporting interval
244  if (k == 0) {
245  continue;
246  }
247 
248  const double t_s = (*m_requested_times)[k - 1];
249  const double t_e = (*m_requested_times)[k];
250 
251  double rate = 0.0;
252  if (N == 1) {
253  // it is the right end-point of the first reporting interval in this time step: count the
254  // contribution from the last time step plus the one from the beginning of this time step
255  const double
256  total_change = m_accumulator + change * (t_e - t0) / (t1 - t0);
257  const double dt = t_e - t_s;
258 
259  rate = total_change / dt;
260 
261  } else {
262  // this reporting interval is completely contained within the time step, so the rate of change
263  // does not depend on its length
264  rate = change / (t1 - t0);
265  }
266 
267  // store computed data in the buffer
268  {
269  m_time.push_back(t_e);
270  m_values.push_back(rate);
271  m_bounds.push_back(t_s);
272  m_bounds.push_back(t_e);
273  }
274 
275  m_accumulator = 0.0;
276  }
277 
278  if (N == 0) {
279  // if this time step contained no requested times we need to add the whole change to the
280  // accumulator
281  m_accumulator += change;
282  } else {
283  // if this time step contained some requested times we need to add the change since the last one
284  // to the accumulator
285  const double dt = t1 - (*m_requested_times)[m_current_time - 1];
286  if (dt > epsilon) {
287  m_accumulator += change * (dt / (t1 - t0));
288  }
289  }
290 }
291 
292 void TSDiagnostic::update(double t0, double t1) {
293  this->update_impl(t0, t1);
294 }
295 
296 void TSSnapshotDiagnostic::update_impl(double t0, double t1) {
297  static const double epsilon = 1e-4; // seconds
298 
299  if (fabs(t1 - t0) < epsilon) {
300  return;
301  }
302 
303  assert(t1 > t0);
304 
305  evaluate(t0, t1, this->compute());
306 }
307 
308 void TSRateDiagnostic::update_impl(double t0, double t1) {
309  const double v = this->compute();
310 
311  if (m_v_previous_set) {
312  assert(t1 > t0);
313  evaluate(t0, t1, v - m_v_previous);
314  }
315 
316  m_v_previous = v;
317  m_v_previous_set = true;
318 }
319 
320 void TSFluxDiagnostic::update_impl(double t0, double t1) {
321  static const double epsilon = 1e-4; // seconds
322 
323  if (fabs(t1 - t0) < epsilon) {
324  return;
325  }
326 
327  assert(t1 > t0);
328 
329  evaluate(t0, t1, this->compute());
330 }
331 
333 
334  if (m_time.empty()) {
335  return;
336  }
337 
338  std::string dimension_name = m_dimension.get_name();
339 
341  io::PISM_READWRITE); // OK to use netcdf3
342 
343  unsigned int len = file.dimension_length(dimension_name);
344 
345  if (len > 0) {
346  double last_time = vector_max(file.read_dimension(dimension_name));
347 
348  if (last_time < m_time.front()) {
349  m_start = len;
350  }
351  }
352 
353  auto time_name = m_config->get_string("time.dimension_name");
354 
355  if (len == m_start) {
357  io::define_time_bounds(m_time_bounds, time_name, "nv", file, io::PISM_DOUBLE);
358 
361  }
362 
365 
366  m_start += m_time.size();
367 
368  {
369  m_time.clear();
370  m_bounds.clear();
371  m_values.clear();
372  }
373 }
374 
375 void TSDiagnostic::init(const File &output_file,
376  std::shared_ptr<std::vector<double>> requested_times) {
377  m_output_filename = output_file.filename();
378 
379  m_requested_times = std::move(requested_times);
380 
381  // Get the number of records in the file (for appending):
382  m_start = output_file.dimension_length(m_dimension.get_name());
383 }
384 
386  return m_variable;
387 }
388 
389 } // end of namespace pism
void write_state(const File &output) const
Definition: Diagnostic.cc:89
virtual void write_state_impl(const File &output) const
Definition: Diagnostic.cc:104
virtual void reset_impl()
Definition: Diagnostic.cc:52
void define(const File &file, io::Type default_type) const
Definition: Diagnostic.cc:120
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:116
double to_internal(double x) const
Definition: Diagnostic.cc:59
virtual std::shared_ptr< array::Array > compute_impl() const =0
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:120
virtual void define_state_impl(const File &output) const
Definition: Diagnostic.cc:99
void init(const File &input, unsigned int time)
Definition: Diagnostic.cc:81
virtual void update_impl(double dt)
Definition: Diagnostic.cc:43
Diagnostic(std::shared_ptr< const Grid > g)
Definition: Diagnostic.cc:31
double to_external(double x) const
Definition: Diagnostic.cc:69
void define_state(const File &output) const
Definition: Diagnostic.cc:85
virtual void define_impl(const File &file, io::Type default_type) const
Define NetCDF variables corresponding to a diagnostic quantity.
Definition: Diagnostic.cc:125
unsigned int n_variables() const
Get the number of NetCDF variables corresponding to a diagnostic quantity.
Definition: Diagnostic.cc:77
void update(double dt)
Definition: Diagnostic.cc:39
std::shared_ptr< const Grid > m_grid
the grid
Definition: Diagnostic.hh:114
virtual void init_impl(const File &input, unsigned int time)
Definition: Diagnostic.cc:93
std::shared_ptr< array::Array > compute() const
Compute a diagnostic quantity and return a pointer to a newly-allocated Array.
Definition: Diagnostic.cc:131
SpatialVariableMetadata & metadata(unsigned int N=0)
Get a metadata object corresponding to variable number N.
Definition: Diagnostic.cc:110
std::string filename() const
Definition: File.cc:307
unsigned int dimension_length(const std::string &name) const
Get the length of a dimension.
Definition: File.cc:454
std::vector< double > read_dimension(const std::string &name) const
Get dimension data (a coordinate variable).
Definition: File.cc:603
High-level PISM I/O class.
Definition: File.hh:56
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
virtual ~TSDiagnostic()
Definition: Diagnostic.cc:167
std::vector< double > m_values
Definition: Diagnostic.hh:327
void set_units(const std::string &units, const std::string &output_units)
Definition: Diagnostic.cc:171
std::shared_ptr< std::vector< double > > m_requested_times
requested times
Definition: Diagnostic.hh:330
VariableMetadata m_dimension
Definition: Diagnostic.hh:321
unsigned int m_current_time
index into m_times
Definition: Diagnostic.hh:332
unsigned int m_start
starting index used when flushing the buffer
Definition: Diagnostic.hh:338
VariableMetadata m_time_bounds
Definition: Diagnostic.hh:322
std::vector< double > m_bounds
Definition: Diagnostic.hh:326
size_t m_buffer_size
size of the buffer used to store data
Definition: Diagnostic.hh:340
VariableMetadata m_variable
Definition: Diagnostic.hh:320
TSDiagnostic(std::shared_ptr< const Grid > g, const std::string &name)
Definition: Diagnostic.cc:145
std::string m_output_filename
Definition: Diagnostic.hh:336
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:313
std::vector< double > m_time
Definition: Diagnostic.hh:325
virtual double compute()=0
void init(const File &output_file, std::shared_ptr< std::vector< double > > requested_times)
Definition: Diagnostic.cc:375
const VariableMetadata & metadata() const
Definition: Diagnostic.cc:385
virtual void update_impl(double t0, double t1)=0
std::shared_ptr< const Grid > m_grid
the grid
Definition: Diagnostic.hh:311
void update(double t0, double t1)
Definition: Diagnostic.cc:292
PISM's scalar time-series diagnostics.
Definition: Diagnostic.hh:278
void update_impl(double t0, double t1)
Definition: Diagnostic.cc:320
TSFluxDiagnostic(std::shared_ptr< const Grid > g, const std::string &name)
Definition: Diagnostic.cc:190
double m_accumulator
accumulator of changes (used to compute rates of change)
Definition: Diagnostic.hh:370
void update_impl(double t0, double t1)
Definition: Diagnostic.cc:308
TSRateDiagnostic(std::shared_ptr< const Grid > g, const std::string &name)
Definition: Diagnostic.cc:185
void evaluate(double t0, double t1, double change)
Definition: Diagnostic.cc:224
double m_v_previous
last two values, used to compute the change during a time step
Definition: Diagnostic.hh:377
Scalar diagnostic reporting the rate of change of a quantity modeled by PISM.
Definition: Diagnostic.hh:364
void evaluate(double t0, double t1, double v)
Definition: Diagnostic.cc:195
TSSnapshotDiagnostic(std::shared_ptr< const Grid > g, const std::string &name)
Definition: Diagnostic.cc:180
void update_impl(double t0, double t1)
Definition: Diagnostic.cc:296
std::string get_name() const
#define PISM_ERROR_LOCATION
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READWRITE
open an existing file for reading and writing
Definition: IO_Flags.hh:74
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
void define_time_bounds(const VariableMetadata &var, const std::string &dimension_name, const std::string &bounds_name, const File &file, io::Type nctype)
Definition: io_helpers.cc:951
void write_time_bounds(const File &file, const VariableMetadata &metadata, size_t t_start, const std::vector< double > &data)
Definition: io_helpers.cc:1057
void define_spatial_variable(const SpatialVariableMetadata &metadata, const Grid &grid, const File &file, io::Type default_type)
Define a NetCDF variable corresponding to a VariableMetadata object.
Definition: io_helpers.cc:467
void write_timeseries(const File &file, const VariableMetadata &metadata, size_t t_start, const std::vector< double > &data)
Write a time-series data to a file.
Definition: io_helpers.cc:926
void define_timeseries(const VariableMetadata &var, const std::string &dimension_name, const File &file, io::Type nctype)
Define a NetCDF variable corresponding to a time-series.
Definition: io_helpers.cc:839
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:70
static const double k
Definition: exactTestP.cc:42
std::string join(const std::vector< std::string > &strings, const std::string &separator)
Concatenate strings, inserting separator between elements.
double vector_max(const std::vector< double > &input)