tSteadyState.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
---
tSteadyState.cc (10880B)
---
1 /* Copyright (C) 2019 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 "EmptyingProblem.hh"
23
24 #include "pism/util/Time.hh" // m_grid->ctx()->time()->current()
25 #include "pism/util/Profiling.hh"
26
27 /* FIXMEs
28 *
29 * - At some later date I might want to support water input rates coming from the model.
30 * In that case I will have to accumulate input at every time step and then use the
31 * average (in time) when it's time to update the model.
32 *
33 * - IceModel::step() updates ice geometry before calling IceModel::hydrology_step(). This
34 * means that a hydrology model has to be able to provide its outputs before the first
35 * update() call. We save subglacial water flux to ensure that this is true for
36 * re-started runs, but in runs started using bootstrapping the first time step will see
37 * "bootstrapped" hydrology outputs (in this case: zero flux).
38 *
39 * - In this context (i.e. computing water flux using piecewise-in-time forcing data) it
40 * makes sense to update the flux at the beginning of an "update interval". In the case
41 * described above (water input coming from the model) we would want to update at the
42 * *end* of an interval.
43 *
44 */
45
46 namespace pism {
47 namespace hydrology {
48
49 void SteadyState::initialization_message() const {
50 m_log->message(2,
51 "* Initializing the \"steady state\" subglacial hydrology model ...\n");
52 }
53
54 SteadyState::SteadyState(IceGrid::ConstPtr grid)
55 : NullTransport(grid) {
56
57 m_time_name = m_config->get_string("time.dimension_name") + "_hydrology_steady";
58 m_t_last = m_grid->ctx()->time()->current();
59 m_update_interval = m_config->get_number("hydrology.steady.flux_update_interval", "seconds");
60 m_t_eps = 1.0;
61 m_bootstrap = false;
62
63 m_emptying_problem.reset(new EmptyingProblem(grid));
64
65 if (m_config->get_flag("hydrology.add_water_input_to_till_storage")) {
66 throw RuntimeError(PISM_ERROR_LOCATION,
67 "'steady' hydrology requires hydrology.add_water_input_to_till_storage == false");
68 }
69
70 if (m_config->get_string("hydrology.surface_input.file").empty()) {
71 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
72 "'steady' hydrology requires hydrology.surface_input.file");
73 }
74 }
75
76 SteadyState::~SteadyState() {
77 // empty
78 }
79
80 void SteadyState::update_impl(double t, double dt, const Inputs& inputs) {
81 NullTransport::update_impl(t, dt, inputs);
82
83 double t_next = m_t_last + max_timestep(m_t_last).value();
84
85 if (t >= t_next or std::abs(t_next - t) < m_t_eps or
86 m_bootstrap) {
87
88 m_log->message(3, " Updating the steady-state subglacial water flux...\n");
89
90 m_grid->ctx()->profiling().begin("steady_emptying");
91
92 m_emptying_problem->update(*inputs.geometry,
93 inputs.no_model_mask,
94 m_surface_input_rate);
95
96 m_grid->ctx()->profiling().end("steady_emptying");
97 m_Q.copy_from(m_emptying_problem->flux());
98
99 m_t_last = t;
100 m_bootstrap = false;
101 }
102 }
103
104 std::map<std::string, Diagnostic::Ptr> SteadyState::diagnostics_impl() const {
105 auto hydro_diagnostics = NullTransport::diagnostics_impl();
106
107 return combine(m_emptying_problem->diagnostics(), hydro_diagnostics);
108 }
109
110 MaxTimestep SteadyState::max_timestep_impl(double t) const {
111
112 // compute the maximum time step coming from the forcing (water input rate)
113 double dt_forcing = 0.0;
114 if (m_time.size() > 0) {
115
116 // the right end point of the last time interval in the forcing file
117 double t_last = m_time_bounds.back();
118
119 double t_next = 0.0;
120 if (t < m_time[0]) {
121 // Allow stepping until the left end point of the first interval.
122 t_next = m_time[0];
123 } else if (t >= t_last) {
124 // Went past the right end point of the last forcing intervals: no time step
125 // restriction from forcing.
126 t_next = t + m_update_interval;
127 } else {
128 // find the index k such that m_time[k] <= t < m_time[k + 1]
129 size_t k = gsl_interp_bsearch(m_time.data(), t, 0, m_time.size());
130
131 assert(m_time[k] <= t);
132 assert(k + 1 == m_time.size() or t < m_time[k + 1]);
133
134 t_next = m_time_bounds[2 * k + 1];
135
136 if (std::abs(t_next - t) < m_t_eps) {
137 // reached t_next; use the right end point of the next interval
138 if (k + 1 < m_time.size()) {
139 t_next = m_time_bounds[2 * (k + 1) + 1];
140 } else {
141 // No time step restriction from forcing. We set dt_forcing to m_update_interval
142 // because dt_interval below will not exceed this, effectively selecting
143 // dt_interval.
144 t_next = t + m_update_interval;
145 }
146 }
147 }
148
149 dt_forcing = t_next - t;
150
151 assert(dt_forcing > 0.0);
152 } else {
153 dt_forcing = m_update_interval;
154 }
155
156 // compute the maximum time step using the update interval
157 double dt_interval = 0.0;
158 {
159 if (t < m_t_last) {
160 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
161 "time %f is less than the previous time %f",
162 t, m_t_last);
163 }
164
165 // Find the smallest time of the form m_t_last + k * m_update_interval that is greater
166 // than t
167 double k = ceil((t - m_t_last) / m_update_interval);
168
169 double
170 t_next = m_t_last + k * m_update_interval;
171
172 dt_interval = t_next - t;
173
174 if (dt_interval < m_t_eps) {
175 dt_interval = m_update_interval;
176 }
177 }
178
179 double dt = std::min(dt_forcing, dt_interval);
180
181 auto dt_null = NullTransport::max_timestep_impl(t);
182 if (dt_null.finite()) {
183 dt = std::min(dt, dt_null.value());
184 }
185
186 return MaxTimestep(dt, "hydrology 'steady'");
187 }
188
189 void SteadyState::define_model_state_impl(const File& output) const {
190 NullTransport::define_model_state_impl(output);
191
192 if (not output.find_variable(m_time_name)) {
193 output.define_variable(m_time_name, PISM_DOUBLE, {});
194
195 output.write_attribute(m_time_name, "long_name",
196 "time of the last update of the steady state subglacial water flux");
197 output.write_attribute(m_time_name, "calendar", m_grid->ctx()->time()->calendar());
198 output.write_attribute(m_time_name, "units", m_grid->ctx()->time()->CF_units_string());
199 }
200
201 m_Q.define(output);
202 }
203
204 void SteadyState::write_model_state_impl(const File& output) const {
205 NullTransport::write_model_state_impl(output);
206
207 output.write_variable(m_time_name, {0}, {1}, &m_t_last);
208 m_Q.write(output);
209 }
210
211 void SteadyState::restart_impl(const File &input_file, int record) {
212 NullTransport::restart_impl(input_file, record);
213
214 init_time(m_config->get_string("hydrology.surface_input.file"));
215
216 // Read m_t_last
217 {
218 if (input_file.find_variable(m_time_name)) {
219 input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
220 } else {
221 m_t_last = m_grid->ctx()->time()->current();
222 }
223 }
224
225 m_Q.read(input_file, record);
226
227 regrid("hydrology 'steady'", m_Q, REGRID_WITHOUT_REGRID_VARS);
228 }
229
230 void SteadyState::bootstrap_impl(const File &input_file,
231 const IceModelVec2S &ice_thickness) {
232 NullTransport::bootstrap_impl(input_file, ice_thickness);
233
234 init_time(m_config->get_string("hydrology.surface_input.file"));
235
236 // Read m_t_last
237 {
238 if (input_file.find_variable(m_time_name)) {
239 input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
240 } else {
241 m_t_last = m_grid->ctx()->time()->current();
242 }
243 }
244
245 // Read water flux
246 if (input_file.find_variable(m_Q.metadata().get_name())) {
247 // Regrid from the input file.
248 m_Q.regrid(input_file, CRITICAL);
249
250 // Allow regridding from a different file.
251 regrid("hydrology 'steady'", m_Q, REGRID_WITHOUT_REGRID_VARS);
252 } else {
253 // No water flux in the input file; try regridding from a different file.
254 auto n = m_Q.state_counter();
255
256 regrid("hydrology 'steady'", m_Q, REGRID_WITHOUT_REGRID_VARS);
257
258 if (n == m_Q.state_counter()) {
259 // Regridding did not fill m_Q: we need to bootstrap during the first update_impl()
260 // call.
261 m_bootstrap = true;
262 }
263 }
264 }
265
266 void SteadyState::init_impl(const IceModelVec2S &W_till,
267 const IceModelVec2S &W,
268 const IceModelVec2S &P) {
269 NullTransport::init_impl(W_till, W, P);
270
271 m_Q.set(0.0);
272
273 m_bootstrap = true;
274 }
275
276
277 /*!
278 * Read time bounds corresponding to the water input rate in the forcing file.
279 *
280 * These times are used to compute the maximum time step the model can take while still
281 * capturing temporal variability of the forcing.
282 */
283 void SteadyState::init_time(const std::string &input_file) {
284
285 std::string variable_name = "water_input_rate";
286
287 File file(m_grid->com, input_file, PISM_GUESS, PISM_READONLY);
288
289 auto time_name = io::time_dimension(m_grid->ctx()->unit_system(),
290 file, variable_name);
291
292 if (time_name.empty()) {
293 // Water input rate is time-independent. m_time and m_time_bounds are left empty.
294 return;
295 }
296
297 std::string bounds_name = file.read_text_attribute(time_name, "bounds");
298
299 if (bounds_name.empty()) {
300 // no time bounds attribute
301 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
302 "Variable '%s' does not have the time_bounds attribute.\n"
303 "Cannot use time-dependent water input rate without time bounds.",
304 time_name.c_str());
305 }
306
307 // read time bounds data from a file
308 TimeBoundsMetadata tb(bounds_name, time_name, m_grid->ctx()->unit_system());
309 tb.set_string("units", m_grid->ctx()->time()->units_string());
310
311 io::read_time_bounds(file, tb, *m_grid->ctx()->time(),
312 *m_log, m_time_bounds);
313
314 // time bounds data overrides the time variable: we make t[j] be the
315 // left end-point of the j-th interval
316 m_time.resize(m_time_bounds.size() / 2);
317 for (unsigned int k = 0; k < m_time.size(); ++k) {
318 m_time[k] = m_time_bounds[2*k];
319 }
320
321 if (not is_increasing(m_time)) {
322 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
323 "time bounds in %s are invalid", input_file.c_str());
324 }
325 }
326
327 } // end of namespace hydrology
328 } // end of namespace pism