PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
SteadyState.cc
Go to the documentation of this file.
1 /* Copyright (C) 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 "SteadyState.hh"
21 
22 #include <gsl/gsl_interp.h> // gsl_interp_bsearch
23 
24 #include "EmptyingProblem.hh"
25 
26 #include "pism/util/Time.hh" // m_grid->ctx()->time()->current()
27 #include "pism/util/Profiling.hh"
28 #include "pism/util/Context.hh"
29 
30 /* FIXMEs
31  *
32  * - At some later date I might want to support water input rates coming from the model.
33  * In that case I will have to accumulate input at every time step and then use the
34  * average (in time) when it's time to update the model.
35  *
36  * - IceModel::step() updates ice geometry before calling IceModel::hydrology_step(). This
37  * means that a hydrology model has to be able to provide its outputs before the first
38  * update() call. We save subglacial water flux to ensure that this is true for
39  * re-started runs, but in runs started using bootstrapping the first time step will see
40  * "bootstrapped" hydrology outputs (in this case: zero flux).
41  *
42  * - In this context (i.e. computing water flux using piecewise-in-time forcing data) it
43  * makes sense to update the flux at the beginning of an "update interval". In the case
44  * described above (water input coming from the model) we would want to update at the
45  * *end* of an interval.
46  *
47  */
48 
49 namespace pism {
50 namespace hydrology {
51 
53  m_log->message(2,
54  "* Initializing the \"steady state\" subglacial hydrology model ...\n");
55 }
56 
58  : NullTransport(grid) {
59 
60  m_time_name = m_config->get_string("time.dimension_name") + "_hydrology_steady";
61  m_t_last = m_grid->ctx()->time()->current();
62  m_update_interval = m_config->get_number("hydrology.steady.flux_update_interval", "seconds");
63  m_t_eps = 1.0;
64  m_bootstrap = false;
65 
67 
68  if (m_config->get_flag("hydrology.add_water_input_to_till_storage")) {
70  "'steady' hydrology requires hydrology.add_water_input_to_till_storage == false");
71  }
72 
73  if (m_config->get_string("hydrology.surface_input.file").empty()) {
75  "'steady' hydrology requires hydrology.surface_input.file");
76  }
77 }
78 
79 void SteadyState::update_impl(double t, double dt, const Inputs& inputs) {
80  NullTransport::update_impl(t, dt, inputs);
81 
82  double t_next = m_t_last + max_timestep(m_t_last).value();
83 
84  if (t >= t_next or std::abs(t_next - t) < m_t_eps or
85  m_bootstrap) {
86 
87  m_log->message(3, " Updating the steady-state subglacial water flux...\n");
88 
89  m_grid->ctx()->profiling().begin("steady_emptying");
90 
91  m_emptying_problem->update(*inputs.geometry,
92  inputs.no_model_mask,
94 
95  m_grid->ctx()->profiling().end("steady_emptying");
97 
98  m_t_last = t;
99  m_bootstrap = false;
100  }
101 }
102 
103 std::map<std::string, Diagnostic::Ptr> SteadyState::diagnostics_impl() const {
104  auto hydro_diagnostics = NullTransport::diagnostics_impl();
105 
106  return combine(m_emptying_problem->diagnostics(), hydro_diagnostics);
107 }
108 
110 
111  // compute the maximum time step coming from the forcing (water input rate)
112  double dt_forcing = 0.0;
113  if (m_time.size() > 0) {
114 
115  // the right end point of the last time interval in the forcing file
116  double t_last = m_time_bounds.back();
117 
118  double t_next = 0.0;
119  if (t < m_time[0]) {
120  // Allow stepping until the left end point of the first interval.
121  t_next = m_time[0];
122  } else if (t >= t_last) {
123  // Went past the right end point of the last forcing intervals: no time step
124  // restriction from forcing.
125  t_next = t + m_update_interval;
126  } else {
127  // find the index k such that m_time[k] <= t < m_time[k + 1]
128  size_t k = gsl_interp_bsearch(m_time.data(), t, 0, m_time.size());
129 
130  assert(m_time[k] <= t);
131  assert(k + 1 == m_time.size() or t < m_time[k + 1]);
132 
133  t_next = m_time_bounds[2 * k + 1];
134 
135  if (std::abs(t_next - t) < m_t_eps) {
136  // reached t_next; use the right end point of the next interval
137  if (k + 1 < m_time.size()) {
138  t_next = m_time_bounds[2 * (k + 1) + 1];
139  } else {
140  // No time step restriction from forcing. We set dt_forcing to m_update_interval
141  // because dt_interval below will not exceed this, effectively selecting
142  // dt_interval.
143  t_next = t + m_update_interval;
144  }
145  }
146  }
147 
148  dt_forcing = t_next - t;
149 
150  assert(dt_forcing > 0.0);
151  } else {
152  dt_forcing = m_update_interval;
153  }
154 
155  // compute the maximum time step using the update interval
156  double dt_interval = 0.0;
157  {
158  if (t < m_t_last) {
160  "time %f is less than the previous time %f",
161  t, m_t_last);
162  }
163 
164  // Find the smallest time of the form m_t_last + k * m_update_interval that is greater
165  // than t
166  double k = ceil((t - m_t_last) / m_update_interval);
167 
168  double
169  t_next = m_t_last + k * m_update_interval;
170 
171  dt_interval = t_next - t;
172 
173  if (dt_interval < m_t_eps) {
174  dt_interval = m_update_interval;
175  }
176  }
177 
178  double dt = std::min(dt_forcing, dt_interval);
179 
180  auto dt_null = NullTransport::max_timestep_impl(t);
181  if (dt_null.finite()) {
182  dt = std::min(dt, dt_null.value());
183  }
184 
185  return MaxTimestep(dt, "hydrology 'steady'");
186 }
187 
188 void SteadyState::define_model_state_impl(const File& output) const {
190 
191  if (not output.find_variable(m_time_name)) {
193 
194  output.write_attribute(m_time_name, "long_name",
195  "time of the last update of the steady state subglacial water flux");
196  output.write_attribute(m_time_name, "calendar", m_grid->ctx()->time()->calendar());
197  output.write_attribute(m_time_name, "units", m_grid->ctx()->time()->units_string());
198  }
199 
200  m_Q.define(output);
201 }
202 
203 void SteadyState::write_model_state_impl(const File& output) const {
205 
206  output.write_variable(m_time_name, {0}, {1}, &m_t_last);
207  m_Q.write(output);
208 }
209 
210 void SteadyState::restart_impl(const File &input_file, int record) {
211  NullTransport::restart_impl(input_file, record);
212 
213  init_time(m_config->get_string("hydrology.surface_input.file"));
214 
215  // Read m_t_last
216  {
217  if (input_file.find_variable(m_time_name)) {
218  input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
219  } else {
220  m_t_last = m_grid->ctx()->time()->current();
221  }
222  }
223 
224  m_Q.read(input_file, record);
225 
226  regrid("hydrology 'steady'", m_Q, REGRID_WITHOUT_REGRID_VARS);
227 }
228 
229 void SteadyState::bootstrap_impl(const File &input_file,
230  const IceModelVec2S &ice_thickness) {
231  NullTransport::bootstrap_impl(input_file, ice_thickness);
232 
233  init_time(m_config->get_string("hydrology.surface_input.file"));
234 
235  // Read m_t_last
236  {
237  if (input_file.find_variable(m_time_name)) {
238  input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
239  } else {
240  m_t_last = m_grid->ctx()->time()->current();
241  }
242  }
243 
244  // Read water flux
245  if (input_file.find_variable(m_Q.metadata().get_name())) {
246  // Regrid from the input file.
247  m_Q.regrid(input_file, CRITICAL);
248 
249  // Allow regridding from a different file.
250  regrid("hydrology 'steady'", m_Q, REGRID_WITHOUT_REGRID_VARS);
251  } else {
252  // No water flux in the input file; try regridding from a different file.
253  auto n = m_Q.state_counter();
254 
255  regrid("hydrology 'steady'", m_Q, REGRID_WITHOUT_REGRID_VARS);
256 
257  if (n == m_Q.state_counter()) {
258  // Regridding did not fill m_Q: we need to bootstrap during the first update_impl()
259  // call.
260  m_bootstrap = true;
261  }
262  }
263 }
264 
266  const IceModelVec2S &W,
267  const IceModelVec2S &P) {
268  NullTransport::init_impl(W_till, W, P);
269 
270  m_Q.set(0.0);
271 
272  m_bootstrap = true;
273 }
274 
275 
276 /*!
277  * Read time bounds corresponding to the water input rate in the forcing file.
278  *
279  * These times are used to compute the maximum time step the model can take while still
280  * capturing temporal variability of the forcing.
281  */
282 void SteadyState::init_time(const std::string &input_file) {
283 
284  std::string variable_name = "water_input_rate";
285 
286  File file(m_grid->com, input_file, PISM_GUESS, PISM_READONLY);
287 
288  auto time_name = io::time_dimension(m_grid->ctx()->unit_system(),
289  file, variable_name);
290 
291  if (time_name.empty()) {
292  // Water input rate is time-independent. m_time and m_time_bounds are left empty.
293  return;
294  }
295 
296  std::string bounds_name = file.read_text_attribute(time_name, "bounds");
297 
298  if (bounds_name.empty()) {
299  // no time bounds attribute
301  "Variable '%s' does not have the time_bounds attribute.\n"
302  "Cannot use time-dependent water input rate without time bounds.",
303  time_name.c_str());
304  }
305 
306  // read time bounds data from a file
307  VariableMetadata tb(bounds_name, m_grid->ctx()->unit_system());
308  tb["units"] = m_grid->ctx()->time()->units_string();
309 
311 
312  // time bounds data overrides the time variable: we make t[j] be the
313  // left end-point of the j-th interval
314  m_time.resize(m_time_bounds.size() / 2);
315  for (unsigned int k = 0; k < m_time.size(); ++k) {
316  m_time[k] = m_time_bounds[2*k];
317  }
318 
319  if (not is_increasing(m_time)) {
321  "time bounds in %s are invalid", input_file.c_str());
322  }
323 }
324 
325 } // end of namespace hydrology
326 } // end of namespace pism
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
@ REGRID_WITHOUT_REGRID_VARS
Definition: Component.hh:131
virtual void regrid(const std::string &module_name, IceModelVec &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:151
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
MaxTimestep max_timestep(double t) const
Reports the maximum time-step the model can take at time t.
Definition: Component.cc:175
void read_variable(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, double *ip) const
Definition: File.cc:740
void define_variable(const std::string &name, IO_Type nctype, const std::vector< std::string > &dims) const
Define a variable.
Definition: File.cc:566
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
Definition: File.cc:364
void write_variable(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, const double *op) const
Definition: File.cc:753
void write_attribute(const std::string &var_name, const std::string &att_name, IO_Type nctype, const std::vector< double > &values) const
Write a multiple-valued double attribute.
Definition: File.cc:631
std::string read_text_attribute(const std::string &var_name, const std::string &att_name) const
Get a text attribute.
Definition: File.cc:686
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
void copy_from(const IceModelVec2< T > &source)
Definition: IceModelVec2.hh:97
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
Definition: iceModelVec.cc:838
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void read(const std::string &filename, unsigned int time)
Definition: iceModelVec.cc:833
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
Definition: iceModelVec.cc:523
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
int state_counter() const
Get the object state counter.
Definition: iceModelVec.cc:124
double value() const
Get the value of the maximum time step.
Definition: MaxTimestep.cc:55
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
std::string get_name() const
IceModelVec2S m_surface_input_rate
Definition: Hydrology.hh:184
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Hydrology.cc:470
virtual std::map< std::string, Diagnostic::Ptr > diagnostics_impl() const
Definition: Hydrology.cc:451
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Hydrology.cc:474
const Geometry * geometry
Definition: Hydrology.hh:42
const IceModelVec2Int * no_model_mask
Definition: Hydrology.hh:40
virtual void restart_impl(const File &input_file, int record)
virtual void init_impl(const IceModelVec2S &W_till, const IceModelVec2S &W, const IceModelVec2S &P)
virtual void bootstrap_impl(const File &input_file, const IceModelVec2S &ice_thickness)
virtual void update_impl(double t, double dt, const Inputs &inputs)
Solves an implicit step of a highly-simplified ODE.
virtual MaxTimestep max_timestep_impl(double t) const
double m_update_interval
Update interval in seconds.
Definition: SteadyState.hh:63
void init_time(const std::string &input_file)
Definition: SteadyState.cc:282
void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: SteadyState.cc:203
double m_t_last
time of the last water flux update
Definition: SteadyState.hh:61
void update_impl(double t, double dt, const Inputs &inputs)
Solves an implicit step of a highly-simplified ODE.
Definition: SteadyState.cc:79
SteadyState(IceGrid::ConstPtr g)
Definition: SteadyState.cc:57
std::vector< double > m_time_bounds
Time bounds corresponding to records in the input file.
Definition: SteadyState.hh:72
MaxTimestep max_timestep_impl(double t) const
Definition: SteadyState.cc:109
std::string m_time_name
Name of the variable used to store the last update time.
Definition: SteadyState.hh:67
void restart_impl(const File &input_file, int record)
Definition: SteadyState.cc:210
bool m_bootstrap
Set to true in bootstrap_impl() if update_impl() has to bootstrap m_Q.
Definition: SteadyState.hh:75
std::map< std::string, Diagnostic::Ptr > diagnostics_impl() const
Definition: SteadyState.cc:103
std::shared_ptr< EmptyingProblem > m_emptying_problem
Definition: SteadyState.hh:58
void bootstrap_impl(const File &input_file, const IceModelVec2S &ice_thickness)
Definition: SteadyState.cc:229
void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: SteadyState.cc:188
std::vector< double > m_time
Times corresponding to records in the input file.
Definition: SteadyState.hh:70
void initialization_message() const
Definition: SteadyState.cc:52
void init_impl(const IceModelVec2S &W_till, const IceModelVec2S &W, const IceModelVec2S &P)
Definition: SteadyState.cc:265
double m_t_eps
Temporal resolution to use when checking whether it's time to update.
Definition: SteadyState.hh:65
#define PISM_ERROR_LOCATION
#define n
Definition: exactTestM.c:37
void read_time_bounds(const File &file, const VariableMetadata &metadata, const Logger &log, std::vector< double > &data)
Definition: io_helpers.cc:1170
std::string time_dimension(units::System::Ptr unit_system, const File &file, const std::string &variable_name)
Definition: io_helpers.cc:1544
bool is_increasing(const std::vector< double > &a)
Checks if a vector of doubles is strictly increasing.
@ PISM_DOUBLE
Definition: IO_Flags.hh:38
static const double k
Definition: exactTestP.cc:45
@ CRITICAL
Definition: IO_Flags.hh:70
@ PISM_GUESS
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
T combine(const T &a, const T &b)