tDistributed.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
---
tDistributed.cc (16018B)
---
1 // Copyright (C) 2012-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 #include <algorithm> // std::min, std::max
20
21 #include "Distributed.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/io/File.hh"
26 #include "pism/util/pism_options.hh"
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/util/IceModelVec2CellType.hh"
29 #include "pism/geometry/Geometry.hh"
30
31 namespace pism {
32 namespace hydrology {
33
34 Distributed::Distributed(IceGrid::ConstPtr g)
35 : Routing(g),
36 m_P(m_grid, "bwp", WITH_GHOSTS, 1),
37 m_Pnew(m_grid, "Pnew_internal", WITHOUT_GHOSTS) {
38
39 // additional variables beyond hydrology::Routing
40 m_P.set_attrs("model_state",
41 "pressure of transportable water in subglacial layer",
42 "Pa", "Pa", "", 0);
43 m_P.metadata().set_number("valid_min", 0.0);
44
45 m_Pnew.set_attrs("internal",
46 "new transportable subglacial water pressure during update",
47 "Pa", "Pa", "", 0);
48 m_Pnew.metadata().set_number("valid_min", 0.0);
49 }
50
51 Distributed::~Distributed() {
52 // empty
53 }
54
55 void Distributed::initialization_message() const {
56 m_log->message(2,
57 "* Initializing the distributed, linked-cavities subglacial hydrology model...\n");
58 }
59
60 void Distributed::restart_impl(const File &input_file, int record) {
61 Routing::restart_impl(input_file, record);
62
63 m_P.read(input_file, record);
64
65 regrid("Hydrology", m_P);
66 }
67
68 void Distributed::bootstrap_impl(const File &input_file,
69 const IceModelVec2S &ice_thickness) {
70 Routing::bootstrap_impl(input_file, ice_thickness);
71
72 double bwp_default = m_config->get_number("bootstrapping.defaults.bwp");
73 m_P.regrid(input_file, OPTIONAL, bwp_default);
74
75 regrid("Hydrology", m_P);
76
77 bool init_P_from_steady = m_config->get_flag("hydrology.distributed.init_p_from_steady");
78
79 if (init_P_from_steady) { // if so, just overwrite -i or -bootstrap value of P=bwp
80 m_log->message(2,
81 " initializing P from P(W) formula which applies in steady state\n");
82
83 compute_overburden_pressure(ice_thickness, m_Pover);
84
85 IceModelVec2S sliding_speed(m_grid, "velbase_mag", WITHOUT_GHOSTS);
86 sliding_speed.set_attrs("internal", "basal sliding speed",
87 "m s-1", "m s-1", "", 0);
88
89 std::string filename = m_config->get_string("hydrology.distributed.sliding_speed_file");
90
91 if (filename.empty()) {
92 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
93 "hydrology.distributed.sliding_speed_file is not set");
94 }
95
96 sliding_speed.regrid(filename, CRITICAL);
97
98 P_from_W_steady(m_W, m_Pover, sliding_speed,
99 m_P);
100 }
101 }
102
103 void Distributed::init_impl(const IceModelVec2S &W_till,
104 const IceModelVec2S &W,
105 const IceModelVec2S &P) {
106 Routing::init_impl(W_till, W, P);
107
108 m_P.copy_from(P);
109 }
110
111 void Distributed::define_model_state_impl(const File &output) const {
112 Routing::define_model_state_impl(output);
113 m_P.define(output);
114 }
115
116 void Distributed::write_model_state_impl(const File &output) const {
117 Routing::write_model_state_impl(output);
118 m_P.write(output);
119 }
120
121 std::map<std::string, TSDiagnostic::Ptr> Distributed::ts_diagnostics_impl() const {
122 std::map<std::string, TSDiagnostic::Ptr> result = {
123 // FIXME: add mass-conservation time-series diagnostics
124 };
125 return result;
126 }
127
128 //! Copies the P state variable which is the modeled water pressure.
129 const IceModelVec2S& Distributed::subglacial_water_pressure() const {
130 return m_P;
131 }
132
133 //! Check bounds on P and fail with message if not satisfied. Optionally, enforces the
134 //! upper bound instead of checking it.
135 /*!
136 The bounds are \f$0 \le P \le P_o\f$ where \f$P_o\f$ is the overburden pressure.
137 */
138 void Distributed::check_P_bounds(IceModelVec2S &P,
139 const IceModelVec2S &P_o,
140 bool enforce_upper) {
141
142 IceModelVec::AccessList list{&P, &P_o};
143
144 ParallelSection loop(m_grid->com);
145 try {
146 for (Points p(*m_grid); p; p.next()) {
147 const int i = p.i(), j = p.j();
148
149 if (P(i,j) < 0.0) {
150 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
151 "negative subglacial water pressure\n"
152 "P(%d, %d) = %.6f Pa",
153 i, j, P(i, j));
154 }
155
156 if (enforce_upper) {
157 P(i,j) = std::min(P(i,j), P_o(i,j));
158 } else if (P(i,j) > P_o(i,j) + 0.001) {
159 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
160 "subglacial water pressure P = %.16f Pa exceeds\n"
161 "overburden pressure Po = %.16f Pa at (i,j)=(%d,%d)",
162 P(i, j), P_o(i, j), i, j);
163 }
164 }
165 } catch (...) {
166 loop.failed();
167 }
168 loop.check();
169
170 }
171
172
173 //! Compute functional relationship P(W) which applies only in steady state.
174 /*!
175 In steady state in this model, water pressure is determined by a balance of
176 cavitation (opening) caused by sliding and creep closure.
177
178 This will be used in initialization when P is otherwise unknown, and
179 in verification and/or reporting. It is not used during time-dependent
180 model runs. To be more complete, \f$P = P(W,P_o,|v_b|)\f$.
181 */
182 void Distributed::P_from_W_steady(const IceModelVec2S &W,
183 const IceModelVec2S &P_overburden,
184 const IceModelVec2S &sliding_speed,
185 IceModelVec2S &result) {
186
187 const double
188 ice_softness = m_config->get_number("flow_law.isothermal_Glen.ice_softness"),
189 creep_closure_coefficient = m_config->get_number("hydrology.creep_closure_coefficient"),
190 cavitation_opening_coefficient = m_config->get_number("hydrology.cavitation_opening_coefficient"),
191 Glen_exponent = m_config->get_number("stress_balance.sia.Glen_exponent"),
192 Wr = m_config->get_number("hydrology.roughness_scale");
193
194 const double CC = cavitation_opening_coefficient / (creep_closure_coefficient * ice_softness);
195
196 IceModelVec::AccessList list{&W, &P_overburden, &sliding_speed, &result};
197
198 for (Points p(*m_grid); p; p.next()) {
199 const int i = p.i(), j = p.j();
200
201 double sb = pow(CC * sliding_speed(i, j), 1.0 / Glen_exponent);
202 if (W(i, j) == 0.0) {
203 // see P(W) formula in steady state; note P(W) is continuous (in steady
204 // state); these facts imply:
205 if (sb > 0.0) {
206 // no water + cavitation = underpressure
207 result(i, j) = 0.0;
208 } else {
209 // no water + no cavitation = creep repressurizes = overburden
210 result(i, j) = P_overburden(i, j);
211 }
212 } else {
213 double Wratio = std::max(0.0, Wr - W(i, j)) / W(i, j);
214 // in cases where steady state is actually possible this will come out positive, but
215 // otherwise we should get underpressure P=0, and that is what it yields
216 result(i, j) = std::max(0.0, P_overburden(i, j) - sb * pow(Wratio, 1.0 / Glen_exponent));
217 }
218 } // end of the loop over grid points
219 }
220
221 double Distributed::max_timestep_P_diff(double phi0, double dt_diff_w) const {
222 return 2.0 * phi0 * dt_diff_w;
223 }
224
225 void Distributed::update_P(double dt,
226 const IceModelVec2CellType &cell_type,
227 const IceModelVec2S &sliding_speed,
228 const IceModelVec2S &surface_input_rate,
229 const IceModelVec2S &basal_melt_rate,
230 const IceModelVec2S &P_overburden,
231 const IceModelVec2S &Wtill,
232 const IceModelVec2S &Wtill_new,
233 const IceModelVec2S &P,
234 const IceModelVec2S &W,
235 const IceModelVec2Stag &Ws,
236 const IceModelVec2Stag &K,
237 const IceModelVec2Stag &Q,
238 IceModelVec2S &P_new) const {
239
240 const double
241 n = m_config->get_number("stress_balance.sia.Glen_exponent"),
242 A = m_config->get_number("flow_law.isothermal_Glen.ice_softness"),
243 c1 = m_config->get_number("hydrology.cavitation_opening_coefficient"),
244 c2 = m_config->get_number("hydrology.creep_closure_coefficient"),
245 Wr = m_config->get_number("hydrology.roughness_scale"),
246 phi0 = m_config->get_number("hydrology.regularizing_porosity");
247
248 // update Pnew from time step
249 const double
250 CC = (m_rg * dt) / phi0,
251 wux = 1.0 / (m_dx * m_dx),
252 wuy = 1.0 / (m_dy * m_dy);
253
254 IceModelVec::AccessList list{&P, &W, &Wtill, &Wtill_new, &sliding_speed, &Ws,
255 &K, &Q, &surface_input_rate, &basal_melt_rate,
256 &cell_type, &P_overburden, &P_new};
257
258 for (Points p(*m_grid); p; p.next()) {
259 const int i = p.i(), j = p.j();
260
261 auto w = W.star(i, j);
262 double P_o = P_overburden(i, j);
263
264 if (cell_type.ice_free_land(i, j)) {
265 P_new(i, j) = 0.0;
266 } else if (cell_type.ocean(i, j)) {
267 P_new(i, j) = P_o;
268 } else if (w.ij <= 0.0) {
269 P_new(i, j) = P_o;
270 } else {
271 auto q = Q.star(i, j);
272 auto k = K.star(i, j);
273 auto ws = Ws.star(i, j);
274
275 double
276 Open = c1 * sliding_speed(i, j) * std::max(0.0, Wr - w.ij),
277 Close = c2 * A * pow(P_o - P(i, j), n) * w.ij;
278
279 // compute the flux divergence the same way as in update_W()
280 const double divadflux = (q.e - q.w) / m_dx + (q.n - q.s) / m_dy;
281 const double
282 De = m_rg * k.e * ws.e,
283 Dw = m_rg * k.w * ws.w,
284 Dn = m_rg * k.n * ws.n,
285 Ds = m_rg * k.s * ws.s;
286
287 double diffW = (wux * (De * (w.e - w.ij) - Dw * (w.ij - w.w)) +
288 wuy * (Dn * (w.n - w.ij) - Ds * (w.ij - w.s)));
289
290 double divflux = -divadflux + diffW;
291
292 // pressure update equation
293 double Wtill_change = Wtill_new(i, j) - Wtill(i, j);
294 double total_input = surface_input_rate(i, j) + basal_melt_rate(i, j);
295 double ZZ = Close - Open + total_input - Wtill_change / dt;
296
297 P_new(i, j) = P(i, j) + CC * (divflux + ZZ);
298
299 // projection to enforce 0 <= P <= P_o
300 P_new(i, j) = clip(P_new(i, j), 0.0, P_o);
301 }
302 } // end of the loop over grid points
303 }
304
305
306 //! Update the model state variables W,P by running the subglacial hydrology model.
307 /*!
308 Runs the hydrology model from time t to time t + dt. Here [t,dt]
309 is generally on the order of months to years. This hydrology model will take its
310 own shorter time steps, perhaps hours to weeks.
311 */
312 void Distributed::update_impl(double t, double dt, const Inputs& inputs) {
313
314 ice_bottom_surface(*inputs.geometry, m_bottom_surface);
315
316 double
317 ht = t,
318 hdt = 0.0;
319
320 const double
321 t_final = t + dt,
322 dt_max = m_config->get_number("hydrology.maximum_time_step", "seconds"),
323 phi0 = m_config->get_number("hydrology.regularizing_porosity");
324
325 m_Qstag_average.set(0.0);
326
327 // make sure W,P have valid ghosts before starting hydrology steps
328 m_W.update_ghosts();
329 m_P.update_ghosts();
330
331 #if (Pism_DEBUG==1)
332 double tillwat_max = m_config->get_number("hydrology.tillwat_max");
333 #endif
334
335 unsigned int step_counter = 0;
336 for (; ht < t_final; ht += hdt) {
337 step_counter++;
338
339 #if (Pism_DEBUG==1)
340 double huge_number = 1e6;
341 check_bounds(m_W, huge_number);
342 check_bounds(m_Wtill, tillwat_max);
343 #endif
344
345 // note that ice dynamics can change overburden pressure, so we can only check P
346 // bounds if thk has not changed; if thk could have just changed, such as in the first
347 // time through the current loop, we enforce them
348 bool enforce_upper = (step_counter == 1);
349 check_P_bounds(m_P, m_Pover, enforce_upper);
350
351 water_thickness_staggered(m_W,
352 inputs.geometry->cell_type,
353 m_Wstag);
354
355 double maxKW = 0.0;
356 compute_conductivity(m_Wstag,
357 subglacial_water_pressure(),
358 m_bottom_surface,
359 m_Kstag, maxKW);
360
361 compute_velocity(m_Wstag,
362 subglacial_water_pressure(),
363 m_bottom_surface,
364 m_Kstag,
365 inputs.no_model_mask,
366 m_Vstag);
367
368 // to get Q, W needs valid ghosts
369 advective_fluxes(m_Vstag, m_W, m_Qstag);
370
371 m_Qstag_average.add(hdt, m_Qstag);
372
373 {
374 const double
375 dt_cfl = max_timestep_W_cfl(),
376 dt_diff_w = max_timestep_W_diff(maxKW),
377 dt_diff_p = max_timestep_P_diff(phi0, dt_diff_w);
378
379 hdt = std::min(t_final - ht, dt_max);
380 hdt = std::min(hdt, dt_cfl);
381 hdt = std::min(hdt, dt_diff_w);
382 hdt = std::min(hdt, dt_diff_p);
383 }
384
385 m_log->message(3, " hydrology step %05d, dt = %f s\n", step_counter, hdt);
386
387 // update Wtillnew from Wtill and input_rate
388 update_Wtill(hdt,
389 m_Wtill,
390 m_surface_input_rate,
391 m_basal_melt_rate,
392 m_Wtillnew);
393 // remove water in ice-free areas and account for changes
394 enforce_bounds(inputs.geometry->cell_type,
395 inputs.no_model_mask,
396 0.0, // do not limit maximum thickness
397 m_Wtillnew,
398 m_grounded_margin_change,
399 m_grounding_line_change,
400 m_conservation_error_change,
401 m_no_model_mask_change);
402
403 update_P(hdt,
404 inputs.geometry->cell_type,
405 *inputs.ice_sliding_speed,
406 m_surface_input_rate,
407 m_basal_melt_rate,
408 m_Pover,
409 m_Wtill, m_Wtillnew,
410 subglacial_water_pressure(),
411 m_W, m_Wstag,
412 m_Kstag, m_Qstag,
413 m_Pnew);
414
415 // update Wnew from W, Wtill, Wtillnew, Wstag, Q, input_rate
416 update_W(hdt,
417 m_surface_input_rate,
418 m_basal_melt_rate,
419 m_W, m_Wstag,
420 m_Wtill, m_Wtillnew,
421 m_Kstag, m_Qstag,
422 m_Wnew);
423 // remove water in ice-free areas and account for changes
424 enforce_bounds(inputs.geometry->cell_type,
425 inputs.no_model_mask,
426 0.0, // do not limit maximum thickness
427 m_Wnew,
428 m_grounded_margin_change,
429 m_grounding_line_change,
430 m_conservation_error_change,
431 m_no_model_mask_change);
432
433 // transfer new into old
434 m_W.copy_from(m_Wnew);
435 m_Wtill.copy_from(m_Wtillnew);
436 m_P.copy_from(m_Pnew);
437 } // end of the time-stepping loop
438
439 staggered_to_regular(inputs.geometry->cell_type, m_Qstag_average,
440 m_config->get_flag("hydrology.routing.include_floating_ice"),
441 m_Q);
442 m_Q.scale(1.0 / dt);
443
444 m_log->message(2,
445 " took %d hydrology sub-steps with average dt = %.6f years (%.6f s)\n",
446 step_counter,
447 units::convert(m_sys, dt/step_counter, "seconds", "years"),
448 dt/step_counter);
449 }
450
451 } // end of namespace hydrology
452 } // end of namespace pism