tLingleClark.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
---
tLingleClark.cc (14676B)
---
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2017, 2018, 2019, 2020 Constantine Khroulev
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 "LingleClark.hh"
20
21 #include "pism/util/io/File.hh"
22 #include "pism/util/Time.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/ConfigInterface.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/Vars.hh"
27 #include "pism/util/MaxTimestep.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/fftw_utilities.hh"
30 #include "LingleClarkSerial.hh"
31
32 namespace pism {
33 namespace bed {
34
35 LingleClark::LingleClark(IceGrid::ConstPtr grid)
36 : BedDef(grid),
37 m_total_displacement(m_grid, "bed_displacement", WITHOUT_GHOSTS),
38 m_relief(m_grid, "bed_relief", WITHOUT_GHOSTS),
39 m_load_thickness(grid, "load_thickness", WITHOUT_GHOSTS),
40 m_elastic_displacement(grid, "elastic_bed_displacement", WITHOUT_GHOSTS) {
41
42 m_time_name = m_config->get_string("time.dimension_name") + "_lingle_clark";
43 m_t_last = m_grid->ctx()->time()->current();
44 m_update_interval = m_config->get_number("bed_deformation.lc.update_interval", "seconds");
45 m_t_eps = 1.0;
46
47 if (m_update_interval < 1.0) {
48 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
49 "invalid bed_deformation.lc.update_interval = %f seconds",
50 m_update_interval);
51 }
52
53 // A work vector. This storage is used to put thickness change on rank 0 and to get the plate
54 // displacement change back.
55 m_total_displacement.set_attrs("internal",
56 "total (viscous and elastic) displacement "
57 "in the Lingle-Clark bed deformation model",
58 "meters", "meters", "", 0);
59
60 m_work0 = m_total_displacement.allocate_proc0_copy();
61
62 m_relief.set_attrs("internal",
63 "bed relief relative to the modeled bed displacement",
64 "meters", "meters", "", 0);
65
66 bool use_elastic_model = m_config->get_flag("bed_deformation.lc.elastic_model");
67
68 m_elastic_displacement.set_attrs("model state",
69 "elastic part of the displacement in the "
70 "Lingle-Clark bed deformation model; "
71 "see :cite:`BLKfastearth`", "meters", "meters", "", 0);
72 m_elastic_displacement0 = m_elastic_displacement.allocate_proc0_copy();
73
74 const int
75 Mx = m_grid->Mx(),
76 My = m_grid->My(),
77 Z = m_config->get_number("bed_deformation.lc.grid_size_factor"),
78 Nx = Z*(Mx - 1) + 1,
79 Ny = Z*(My - 1) + 1;
80
81 const double
82 Lx = Z * (m_grid->x0() - m_grid->x(0)),
83 Ly = Z * (m_grid->y0() - m_grid->y(0));
84
85 m_extended_grid = IceGrid::Shallow(m_grid->ctx(),
86 Lx, Ly,
87 m_grid->x0(), m_grid->y0(),
88 Nx, Ny, CELL_CORNER, NOT_PERIODIC);
89
90 m_viscous_displacement.create(m_extended_grid,
91 "viscous_bed_displacement", WITHOUT_GHOSTS);
92 m_viscous_displacement.set_attrs("model state",
93 "bed displacement in the viscous half-space "
94 "bed deformation model; "
95 "see BuelerLingleBrown",
96 "meters", "meters", "", 0);
97
98 // coordinate variables of the extended grid should have different names
99 m_viscous_displacement.metadata().get_x().set_name("x_lc");
100 m_viscous_displacement.metadata().get_y().set_name("y_lc");
101
102 // do not point to auxiliary coordinates "lon" and "lat".
103 m_viscous_displacement.metadata().set_string("coordinates", "");
104
105 m_viscous_displacement0 = m_viscous_displacement.allocate_proc0_copy();
106
107 ParallelSection rank0(m_grid->com);
108 try {
109 if (m_grid->rank() == 0) {
110 m_serial_model.reset(new LingleClarkSerial(m_log, *m_config, use_elastic_model,
111 Mx, My,
112 m_grid->dx(), m_grid->dy(),
113 Nx, Ny));
114 }
115 } catch (...) {
116 rank0.failed();
117 }
118 rank0.check();
119 }
120
121 LingleClark::~LingleClark() {
122 // empty
123 }
124
125 /*!
126 * Initialize the model by computing the viscous bed displacement using uplift and the elastic
127 * response using ice thickness.
128 *
129 * Then compute the bed relief as the difference between bed elevation and total bed displacement.
130 *
131 * This method has to initialize m_viscous_displacement, m_elastic_displacement,
132 * m_total_displacement, and m_relief.
133 */
134 void LingleClark::bootstrap_impl(const IceModelVec2S &bed_elevation,
135 const IceModelVec2S &bed_uplift,
136 const IceModelVec2S &ice_thickness,
137 const IceModelVec2S &sea_level_elevation) {
138 m_t_last = m_grid->ctx()->time()->current();
139
140 m_topg_last.copy_from(bed_elevation);
141
142 compute_load(bed_elevation, ice_thickness, sea_level_elevation,
143 m_load_thickness);
144
145 petsc::Vec::Ptr thickness0 = m_load_thickness.allocate_proc0_copy();
146
147 // initialize the plate displacement
148 {
149 bed_uplift.put_on_proc0(*m_work0);
150 m_load_thickness.put_on_proc0(*thickness0);
151
152 ParallelSection rank0(m_grid->com);
153 try {
154 if (m_grid->rank() == 0) {
155 PetscErrorCode ierr = 0;
156
157 m_serial_model->bootstrap(*thickness0, *m_work0);
158
159 ierr = VecCopy(m_serial_model->total_displacement(), *m_work0);
160 PISM_CHK(ierr, "VecCopy");
161
162 ierr = VecCopy(m_serial_model->viscous_displacement(), *m_viscous_displacement0);
163 PISM_CHK(ierr, "VecCopy");
164
165 ierr = VecCopy(m_serial_model->elastic_displacement(), *m_elastic_displacement0);
166 PISM_CHK(ierr, "VecCopy");
167 }
168 } catch (...) {
169 rank0.failed();
170 }
171 rank0.check();
172 }
173
174 m_viscous_displacement.get_from_proc0(*m_viscous_displacement0);
175
176 m_elastic_displacement.get_from_proc0(*m_elastic_displacement0);
177
178 m_total_displacement.get_from_proc0(*m_work0);
179
180 // compute bed relief
181 m_topg.add(-1.0, m_total_displacement, m_relief);
182 }
183
184 /*!
185 * Return the load response matrix for the elastic response.
186 *
187 * This method is used for testing only.
188 */
189 IceModelVec2S::Ptr LingleClark::elastic_load_response_matrix() const {
190 IceModelVec2S::Ptr result(new IceModelVec2S(m_extended_grid, "lrm", WITHOUT_GHOSTS));
191
192 int
193 Nx = m_extended_grid->Mx(),
194 Ny = m_extended_grid->My();
195
196 auto lrm0 = result->allocate_proc0_copy();
197
198 {
199 ParallelSection rank0(m_grid->com);
200 try {
201 if (m_grid->rank() == 0) {
202 std::vector<std::complex<double> > array(Nx * Ny);
203
204 m_serial_model->compute_load_response_matrix((fftw_complex*)array.data());
205
206 get_real_part((fftw_complex*)array.data(), 1.0, Nx, Ny, Nx, Ny, 0, 0, *lrm0);
207 }
208 } catch (...) {
209 rank0.failed();
210 }
211 rank0.check();
212 }
213
214 result->get_from_proc0(*lrm0);
215
216 return result;
217 }
218
219 /*! Initialize the Lingle-Clark bed deformation model.
220 *
221 * Inputs:
222 *
223 * - bed topography,
224 * - ice thickness,
225 * - plate displacement (either read from a file or bootstrapped using uplift) and
226 * possibly re-gridded.
227 */
228 void LingleClark::init_impl(const InputOptions &opts, const IceModelVec2S &ice_thickness,
229 const IceModelVec2S &sea_level_elevation) {
230 m_log->message(2, "* Initializing the Lingle-Clark bed deformation model...\n");
231
232 if (opts.type == INIT_RESTART or opts.type == INIT_BOOTSTRAP) {
233 File input_file(m_grid->com, opts.filename, PISM_NETCDF3, PISM_READONLY);
234
235 if (input_file.find_variable(m_time_name)) {
236 input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
237 } else {
238 m_t_last = m_grid->ctx()->time()->current();
239 }
240 } else {
241 m_t_last = m_grid->ctx()->time()->current();
242 }
243
244 // Initialize bed topography and uplift maps.
245 BedDef::init_impl(opts, ice_thickness, sea_level_elevation);
246
247 m_topg_last.copy_from(m_topg);
248
249 if (opts.type == INIT_RESTART) {
250 // Set viscous displacement by reading from the input file.
251 m_viscous_displacement.read(opts.filename, opts.record);
252 // Set elastic displacement by reading from the input file.
253 m_elastic_displacement.read(opts.filename, opts.record);
254 } else if (opts.type == INIT_BOOTSTRAP) {
255 this->bootstrap(m_topg, m_uplift, ice_thickness, sea_level_elevation);
256 } else {
257 // do nothing
258 }
259
260 // Try re-gridding plate_displacement.
261 regrid("Lingle-Clark bed deformation model",
262 m_viscous_displacement, REGRID_WITHOUT_REGRID_VARS);
263 regrid("Lingle-Clark bed deformation model",
264 m_elastic_displacement, REGRID_WITHOUT_REGRID_VARS);
265
266 compute_load(m_topg, ice_thickness, sea_level_elevation,
267 m_load_thickness);
268
269 // Now that viscous displacement and elastic displacement are finally initialized,
270 // put them on rank 0 and initialize the serial model itself.
271 {
272 m_viscous_displacement.put_on_proc0(*m_viscous_displacement0);
273 m_elastic_displacement.put_on_proc0(*m_work0);
274
275 ParallelSection rank0(m_grid->com);
276 try {
277 if (m_grid->rank() == 0) { // only processor zero does the work
278 PetscErrorCode ierr = 0;
279
280 m_serial_model->init(*m_viscous_displacement0, *m_work0);
281
282 ierr = VecCopy(m_serial_model->total_displacement(), *m_work0);
283 PISM_CHK(ierr, "VecCopy");
284 }
285 } catch (...) {
286 rank0.failed();
287 }
288 rank0.check();
289 }
290
291 m_total_displacement.get_from_proc0(*m_work0);
292
293 // compute bed relief
294 m_topg.add(-1.0, m_total_displacement, m_relief);
295 }
296
297 MaxTimestep LingleClark::max_timestep_impl(double t) const {
298
299 if (t < m_t_last) {
300 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
301 "time %f is less than the previous time %f",
302 t, m_t_last);
303 }
304
305 // Find the smallest time of the form m_t_last + k * m_update_interval that is greater
306 // than t
307 double k = ceil((t - m_t_last) / m_update_interval);
308
309 double
310 t_next = m_t_last + k * m_update_interval,
311 dt_max = t_next - t;
312
313 if (dt_max < m_t_eps) {
314 dt_max = m_update_interval;
315 }
316
317 return MaxTimestep(dt_max, "bed_def lc");
318 }
319
320 /*!
321 * Get total bed displacement on the PISM grid.
322 */
323 const IceModelVec2S& LingleClark::total_displacement() const {
324 return m_total_displacement;
325 }
326
327 const IceModelVec2S& LingleClark::viscous_displacement() const {
328 return m_viscous_displacement;
329 }
330
331 const IceModelVec2S& LingleClark::elastic_displacement() const {
332 return m_elastic_displacement;
333 }
334
335 const IceModelVec2S& LingleClark::relief() const {
336 return m_relief;
337 }
338
339 void LingleClark::step(const IceModelVec2S &ice_thickness,
340 const IceModelVec2S &sea_level_elevation,
341 double dt) {
342
343 compute_load(m_topg, ice_thickness, sea_level_elevation,
344 m_load_thickness);
345
346 m_load_thickness.put_on_proc0(*m_work0);
347
348 ParallelSection rank0(m_grid->com);
349 try {
350 if (m_grid->rank() == 0) { // only processor zero does the step
351 PetscErrorCode ierr = 0;
352
353 m_serial_model->step(dt, *m_work0);
354
355 ierr = VecCopy(m_serial_model->total_displacement(), *m_work0);
356 PISM_CHK(ierr, "VecCopy");
357
358 ierr = VecCopy(m_serial_model->viscous_displacement(), *m_viscous_displacement0);
359 PISM_CHK(ierr, "VecCopy");
360
361 ierr = VecCopy(m_serial_model->elastic_displacement(), *m_elastic_displacement0);
362 PISM_CHK(ierr, "VecCopy");
363 }
364 } catch (...) {
365 rank0.failed();
366 }
367 rank0.check();
368
369 m_viscous_displacement.get_from_proc0(*m_viscous_displacement0);
370
371 m_elastic_displacement.get_from_proc0(*m_elastic_displacement0);
372
373 m_total_displacement.get_from_proc0(*m_work0);
374
375 // Update bed elevation using bed displacement and relief.
376 {
377 m_total_displacement.add(1.0, m_relief, m_topg);
378 // Increment the topg state counter. SIAFD relies on this!
379 m_topg.inc_state_counter();
380 }
381
382 //! Finally, we need to update bed uplift and topg_last.
383 compute_uplift(m_topg, m_topg_last, dt, m_uplift);
384 m_topg_last.copy_from(m_topg);
385 }
386
387 //! Update the Lingle-Clark bed deformation model.
388 void LingleClark::update_impl(const IceModelVec2S &ice_thickness,
389 const IceModelVec2S &sea_level_elevation,
390 double t, double dt) {
391
392 double
393 t_next = m_t_last + m_update_interval,
394 t_final = t + dt;
395
396 if (t_final < m_t_last) {
397 throw RuntimeError(PISM_ERROR_LOCATION, "cannot go back in time");
398 }
399
400 if (std::abs(t_next - t_final) < m_t_eps) { // reached the next update time
401 double dt_beddef = t_final - m_t_last;
402 step(ice_thickness, sea_level_elevation, dt_beddef);
403 m_t_last = t_final;
404 }
405 }
406
407 void LingleClark::define_model_state_impl(const File &output) const {
408 BedDef::define_model_state_impl(output);
409 m_viscous_displacement.define(output);
410 m_elastic_displacement.define(output);
411
412 if (not output.find_variable(m_time_name)) {
413 output.define_variable(m_time_name, PISM_DOUBLE, {});
414
415 output.write_attribute(m_time_name, "long_name",
416 "time of the last update of the Lingle-Clark bed deformation model");
417 output.write_attribute(m_time_name, "calendar", m_grid->ctx()->time()->calendar());
418 output.write_attribute(m_time_name, "units", m_grid->ctx()->time()->CF_units_string());
419 }
420 }
421
422 void LingleClark::write_model_state_impl(const File &output) const {
423 BedDef::write_model_state_impl(output);
424
425 m_viscous_displacement.write(output);
426 m_elastic_displacement.write(output);
427
428 output.write_variable(m_time_name, {0}, {1}, &m_t_last);
429 }
430
431 DiagnosticList LingleClark::diagnostics_impl() const {
432 DiagnosticList result = {
433 {"viscous_bed_displacement", Diagnostic::wrap(m_viscous_displacement)},
434 {"elastic_bed_displacement", Diagnostic::wrap(m_elastic_displacement)},
435 };
436 return combine(result, BedDef::diagnostics_impl());
437 }
438
439 } // end of namespace bed
440 } // end of namespace pism