tpismi.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tpismi.py (23418B)
---
1 #! /usr/bin/env python3
2 #
3 # Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 David Maxwell and Constantine Khroulev
4 #
5 # This file is part of PISM.
6 #
7 # PISM is free software; you can redistribute it and/or modify it under the
8 # terms of the GNU General Public License as published by the Free Software
9 # Foundation; either version 3 of the License, or (at your option) any later
10 # version.
11 #
12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 # details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with PISM; if not, write to the Free Software
19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21 # try to start coverage
22 try: # pragma: no cover
23 import coverage
24 cov = coverage.coverage(branch=True)
25 try:
26 # try to load coverage data and ignore failures
27 cov.load()
28 except:
29 pass
30 cov.start()
31 except ImportError: # pragma: no cover
32 pass
33
34 import PISM
35 import PISM.invert.ssa
36 from PISM.logging import logMessage
37 from PISM.util import convert
38 import numpy as np
39 import sys, os, math
40
41
42 class SSAForwardRun(PISM.invert.ssa.SSAForwardRunFromInputFile):
43
44 def write(self, filename, append=False):
45 if not append:
46 PISM.invert.ssa.SSAForwardRunFromInputFile.write(self, filename)
47 else:
48 grid = self.grid
49 vecs = self.modeldata.vecs
50
51 pio = PISM.File(grid.com, filename, PISM.PISM_NETCDF3, PISM.PISM_READWRITE)
52
53 self.modeldata.vecs.write(filename)
54 pio.close()
55
56
57 class InvSSAPlotListener(PISM.invert.listener.PlotListener):
58
59 def __init__(self, grid, Vmax):
60 PISM.invert.listener.PlotListener.__init__(self, grid)
61 self.Vmax = Vmax
62 self.l2_weight = None
63 self.l2_weight_init = False
64
65 def __call__(self, inverse_solver, count, data):
66
67 if not self.l2_weight_init:
68 vecs = inverse_solver.ssarun.modeldata.vecs
69 if vecs.has('vel_misfit_weight'):
70 self.l2_weight = self.toproczero(vecs.vel_misfit_weight)
71 self.l2_weight_init = True
72
73 method = inverse_solver.method
74
75 r = self.toproczero(data.residual)
76 Td = None
77 if 'T_zeta_step' in data:
78 Td = self.toproczero(data.T_zeta_step)
79 TStarR = None
80 if 'TStar_residual' in data:
81 TStarR = self.toproczero(data.TStar_residual)
82 d = None
83 if 'zeta_step' in data:
84 d = self.toproczero(data.zeta_step)
85 zeta = self.toproczero(data.zeta)
86
87 secpera = convert(1.0, "year", "second")
88
89 if self.grid.rank() == 0:
90 import matplotlib.pyplot as pp
91
92 pp.figure(self.figure())
93
94 l2_weight = self.l2_weight
95
96 pp.clf()
97
98 V = self.Vmax
99
100 pp.subplot(2, 3, 1)
101 if l2_weight is not None:
102 rx = l2_weight * r[0, :, :] * secpera
103 else:
104 rx = r[0, :, :] * secpera
105 rx = np.maximum(rx, -V)
106 rx = np.minimum(rx, V)
107 pp.imshow(rx, origin='lower', interpolation='nearest')
108 pp.colorbar()
109 pp.title('r_x')
110 pp.jet()
111
112 pp.subplot(2, 3, 4)
113 if l2_weight is not None:
114 ry = l2_weight * r[1, :, :] * secpera
115 else:
116 ry = r[1, :, :] * secpera
117 ry = np.maximum(ry, -V)
118 ry = np.minimum(ry, V)
119 pp.imshow(ry, origin='lower', interpolation='nearest')
120 pp.colorbar()
121 pp.title('r_y')
122 pp.jet()
123
124 if method == 'ign':
125 pp.subplot(2, 3, 2)
126 Tdx = Td[0, :, :] * secpera
127 pp.imshow(Tdx, origin='lower', interpolation='nearest')
128 pp.colorbar()
129 pp.title('Td_x')
130 pp.jet()
131
132 pp.subplot(2, 3, 5)
133 Tdy = Td[1, :, :] * secpera
134 pp.imshow(Tdy, origin='lower', interpolation='nearest')
135 pp.colorbar()
136 pp.title('Td_y')
137 pp.jet()
138 elif method == 'sd' or method == 'nlcg':
139 pp.subplot(2, 3, 2)
140 pp.imshow(TStarR, origin='lower', interpolation='nearest')
141 pp.colorbar()
142 pp.title('TStarR')
143 pp.jet()
144
145 if d is not None:
146 d *= -1
147 pp.subplot(2, 3, 3)
148 pp.imshow(d, origin='lower', interpolation='nearest')
149
150 # colorbar does a divide by zero if 'd' is all zero,
151 # as it will be at the start of iteration zero.
152 # The warning message is a distraction, so we suppress it.
153 import warnings
154 with warnings.catch_warnings():
155 warnings.simplefilter("ignore")
156 pp.colorbar()
157 pp.jet()
158
159 pp.title('-zeta_step')
160
161 pp.subplot(2, 3, 6)
162 pp.imshow(zeta, origin='lower', interpolation='nearest')
163 pp.colorbar()
164 pp.jet()
165 pp.title('zeta')
166
167 pp.ion()
168 pp.draw()
169 pp.show()
170
171
172 class InvSSALinPlotListener(PISM.invert.listener.PlotListener):
173
174 def __init__(self, grid, Vmax):
175 PISM.invert.listener.PlotListener.__init__(self, grid)
176 self.Vmax = Vmax
177 self.l2_weight = None
178 self.l2_weight_init = False
179
180 def __call__(self, inverse_solver, count, data):
181 # On the first go-around, extract the l2_weight vector onto
182 # processor zero.
183 if self.l2_weight_init == False:
184 vecs = inverse_solver.ssarun.modeldata.vecs
185 self.l2_weight = self.toproczero(vecs.vel_misfit_weight)
186 self.l2_init = True
187
188 l2_weight = self.l2_weight
189 r = self.toproczero(data.r)
190 d = self.toproczero(data.d)
191
192 if self.grid.rank() == 0:
193 import matplotlib.pyplot as pp
194 pp.figure(self.figure())
195 pp.clf()
196
197 V = self.Vmax
198 pp.subplot(1, 3, 1)
199 rx = l2_weight * r[0, :, :]
200 rx = np.maximum(rx, -V)
201 rx = np.minimum(rx, V)
202 pp.imshow(rx, origin='lower', interpolation='nearest')
203 pp.colorbar()
204 pp.title('ru')
205 pp.jet()
206
207 pp.subplot(1, 3, 2)
208 ry = l2_weight * r[1, :, :]
209 ry = np.maximum(ry, -V)
210 ry = np.minimum(ry, V)
211 pp.imshow(ry, origin='lower', interpolation='nearest')
212 pp.colorbar()
213 pp.title('rv')
214 pp.jet()
215
216 d *= -1
217 pp.subplot(1, 3, 3)
218 pp.imshow(d, origin='lower', interpolation='nearest')
219 pp.colorbar()
220 pp.jet()
221 pp.title('-d')
222
223 pp.ion()
224 pp.show()
225
226
227 def adjustTauc(mask, tauc):
228 """Where ice is floating or land is ice-free, tauc should be adjusted to have some preset default values."""
229
230 logMessage(" Adjusting initial estimate of 'tauc' to match PISM model for floating ice and ice-free bedrock.\n")
231
232 grid = mask.grid()
233 high_tauc = grid.ctx().config().get_number("basal_yield_stress.ice_free_bedrock")
234
235 with PISM.vec.Access(comm=tauc, nocomm=mask):
236 for (i, j) in grid.points():
237 if mask.ocean(i, j):
238 tauc[i, j] = 0
239 elif mask.ice_free(i, j):
240 tauc[i, j] = high_tauc
241
242
243 def createDesignVec(grid, design_var, name=None, **kwargs):
244 if name is None:
245 name = design_var
246 if design_var == "tauc":
247 design_vec = PISM.model.createYieldStressVec(grid, name=name, **kwargs)
248 elif design_var == "hardav":
249 design_vec = PISM.model.createAveragedHardnessVec(grid, name=name, **kwargs)
250 else:
251 raise ValueError("Unknown design variable %s" % design_var)
252 return design_vec
253
254
255 # Main code starts here
256 def run():
257 context = PISM.Context()
258 config = context.config
259 com = context.com
260 PISM.set_abort_on_sigint(True)
261
262 WIDE_STENCIL = int(config.get_number("grid.max_stencil_width"))
263
264 usage = \
265 """ pismi.py [-i IN.nc [-o OUT.nc]]/[-a INOUT.nc] [-inv_data inv_data.nc] [-inv_forward model]
266 [-inv_design design_var] [-inv_method meth]
267 where:
268 -i IN.nc is input file in NetCDF format: contains PISM-written model state
269 -o OUT.nc is output file in NetCDF format to be overwritten
270 -a INOUT.nc is input/output file in NetCDF format to be appended to
271 -inv_data inv_data.nc is data file containing extra inversion data (e.g. observed surface velocities)
272 -inv_forward model forward model: only 'ssa' supported
273 -inv_design design_var design variable name; one of 'tauc'/'hardav' for SSA inversions
274 -inv_method meth algorithm for inversion [sd,nlcg,ign,tikhonov_lmvm]
275
276 notes:
277 * only one of -i/-a is allowed; both specify the input file
278 * only one of -o/-a is allowed; both specify the output file
279 * if -o is used, only the variables involved in inversion are written to the output file.
280 * if -a is used, the varaibles involved in inversion are appended to the given file. No
281 original variables in the file are changed.
282 """
283
284 append_mode = False
285
286 input_filename = config.get_string("input.file")
287 if len(input_filename) == 0:
288 input_filename = None
289
290 append = PISM.OptionString("-a", "append file")
291 append_filename = append.value() if append.is_set() else None
292
293 output_filename = config.get_string("output.file_name")
294 if len(output_filename) == 0:
295 output_filename = None
296
297 if (input_filename is None) and (append_filename is None):
298 PISM.verbPrintf(1, com, "\nError: No input file specified. Use one of -i [file.nc] or -a [file.nc].\n")
299 sys.exit(0)
300
301 if (input_filename is not None) and (append_filename is not None):
302 PISM.verbPrintf(1, com, "\nError: Only one of -i/-a is allowed.\n")
303 sys.exit(0)
304
305 if (output_filename is not None) and (append_filename is not None):
306 PISM.verbPrintf(1, com, "\nError: Only one of -a/-o is allowed.\n")
307 sys.edit(0)
308
309 if append_filename is not None:
310 input_filename = append_filename
311 output_filename = append_filename
312 append_mode = True
313
314 inv_data_filename = PISM.OptionString("-inv_data", "inverse data file", input_filename).value()
315
316 do_plotting = PISM.OptionBool("-inv_plot", "perform visualization during the computation")
317 do_final_plot = PISM.OptionBool("-inv_final_plot",
318 "perform visualization at the end of the computation")
319 Vmax = PISM.OptionReal("-inv_plot_vmax", "maximum velocity for plotting residuals", 30)
320
321 design_var = PISM.OptionKeyword("-inv_ssa",
322 "design variable for inversion",
323 "tauc,hardav", "tauc").value()
324 do_pause = PISM.OptionBool("-inv_pause", "pause each iteration")
325
326 do_restart = PISM.OptionBool("-inv_restart", "Restart a stopped computation.")
327 use_design_prior = config.get_flag("inverse.use_design_prior")
328
329 prep_module = PISM.OptionString("-inv_prep_module",
330 "Python module used to do final setup of inverse solver")
331 prep_module = prep_module.value() if prep_module.is_set() else None
332
333 is_regional = PISM.OptionBool("-regional", "Compute SIA/SSA using regional model semantics")
334
335 using_zeta_fixed_mask = config.get_flag("inverse.use_zeta_fixed_mask")
336
337 inv_method = config.get_string("inverse.ssa.method")
338
339 if output_filename is None:
340 output_filename = "pismi_" + os.path.basename(input_filename)
341
342 saving_inv_data = (inv_data_filename != output_filename)
343
344 forward_run = SSAForwardRun(input_filename, inv_data_filename, design_var)
345 forward_run.setup()
346 design_param = forward_run.designVariableParameterization()
347 solver = PISM.invert.ssa.createInvSSASolver(forward_run)
348
349 modeldata = forward_run.modeldata
350 vecs = modeldata.vecs
351 grid = modeldata.grid
352
353 # Determine the prior guess for tauc/hardav. This can be one of
354 # a) tauc/hardav from the input file (default)
355 # b) tauc/hardav_prior from the inv_datafile if -inv_use_design_prior is set
356 design_prior = createDesignVec(grid, design_var, '%s_prior' % design_var)
357 long_name = design_prior.metadata().get_string("long_name")
358 units = design_prior.metadata().get_string("units")
359 design_prior.set_attrs("", "best prior estimate for %s (used for inversion)" % long_name,
360 units, units, "", 0)
361 if PISM.util.fileHasVariable(inv_data_filename, "%s_prior" % design_var) and use_design_prior:
362 PISM.logging.logMessage(" Reading '%s_prior' from inverse data file %s.\n" % (design_var, inv_data_filename))
363 design_prior.regrid(inv_data_filename, critical=True)
364 vecs.add(design_prior, writing=saving_inv_data)
365 else:
366 if not PISM.util.fileHasVariable(input_filename, design_var):
367 PISM.verbPrintf(1, com, "Initial guess for design variable is not available as '%s' in %s.\nYou can provide an initial guess in the inverse data file.\n" % (
368 design_var, input_filename))
369 exit(1)
370 PISM.logging.logMessage("Reading '%s_prior' from '%s' in input file.\n" % (design_var, design_var))
371 design = createDesignVec(grid, design_var)
372 design.regrid(input_filename, True)
373 design_prior.copy_from(design)
374 vecs.add(design_prior, writing=True)
375
376 if using_zeta_fixed_mask:
377 if PISM.util.fileHasVariable(inv_data_filename, "zeta_fixed_mask"):
378 zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
379 zeta_fixed_mask.regrid(inv_data_filename)
380 vecs.add(zeta_fixed_mask)
381 else:
382 if design_var == 'tauc':
383 logMessage(
384 " Computing 'zeta_fixed_mask' (i.e. locations where design variable '%s' has a fixed value).\n" % design_var)
385 zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
386 zeta_fixed_mask.set(1)
387 mask = vecs.mask
388 with PISM.vec.Access(comm=zeta_fixed_mask, nocomm=mask):
389 for (i, j) in grid.points():
390 if mask.grounded_ice(i, j):
391 zeta_fixed_mask[i, j] = 0
392 vecs.add(zeta_fixed_mask)
393
394 adjustTauc(vecs.mask, design_prior)
395 elif design_var == 'hardav':
396 PISM.logging.logPrattle(
397 "Skipping 'zeta_fixed_mask' for design variable 'hardav'; no natural locations to fix its value.")
398 pass
399 else:
400 raise NotImplementedError("Unable to build 'zeta_fixed_mask' for design variable %s.", design_var)
401
402 # Convert design_prior -> zeta_prior
403 zeta_prior = PISM.IceModelVec2S()
404 zeta_prior.create(grid, "zeta_prior", PISM.WITH_GHOSTS, WIDE_STENCIL)
405 design_param.convertFromDesignVariable(design_prior, zeta_prior)
406 vecs.add(zeta_prior, writing=True)
407
408 # Determine the initial guess for zeta. If we are restarting, load it from
409 # the output file. Otherwise, if 'zeta_inv' is in the inverse data file, use it.
410 # If none of the above, copy from 'zeta_prior'.
411 zeta = PISM.IceModelVec2S()
412 zeta.create(grid, "zeta_inv", PISM.WITH_GHOSTS, WIDE_STENCIL)
413 zeta.set_attrs("diagnostic", "zeta_inv", "1", "1", "zeta_inv", 0)
414 if do_restart:
415 # Just to be sure, verify that we have a 'zeta_inv' in the output file.
416 if not PISM.util.fileHasVariable(output_filename, 'zeta_inv'):
417 PISM.verbPrintf(
418 1, com, "Unable to restart computation: file %s is missing variable 'zeta_inv'", output_filename)
419 exit(1)
420 PISM.logging.logMessage(" Inversion starting from 'zeta_inv' found in %s\n" % output_filename)
421 zeta.regrid(output_filename, True)
422
423 elif PISM.util.fileHasVariable(inv_data_filename, 'zeta_inv'):
424 PISM.logging.logMessage(" Inversion starting from 'zeta_inv' found in %s\n" % inv_data_filename)
425 zeta.regrid(inv_data_filename, True)
426 else:
427 zeta.copy_from(zeta_prior)
428
429 vel_ssa_observed = None
430 vel_ssa_observed = PISM.model.create2dVelocityVec(grid, '_ssa_observed', stencil_width=2)
431 if PISM.util.fileHasVariable(inv_data_filename, "u_ssa_observed"):
432 vel_ssa_observed.regrid(inv_data_filename, True)
433 vecs.add(vel_ssa_observed, writing=saving_inv_data)
434 else:
435 if not PISM.util.fileHasVariable(inv_data_filename, "u_surface_observed"):
436 PISM.verbPrintf(
437 1, context.com, "Neither u/v_ssa_observed nor u/v_surface_observed is available in %s.\nAt least one must be specified.\n" % inv_data_filename)
438 exit(1)
439 vel_surface_observed = PISM.model.create2dVelocityVec(grid, '_surface_observed', stencil_width=2)
440 vel_surface_observed.regrid(inv_data_filename, True)
441 vecs.add(vel_surface_observed, writing=saving_inv_data)
442
443 sia_solver = PISM.SIAFD
444 if is_regional:
445 sia_solver = PISM.SIAFD_Regional
446 vel_sia_observed = PISM.sia.computeSIASurfaceVelocities(modeldata, sia_solver)
447
448 vel_sia_observed.metadata(0).set_name('u_sia_observed')
449 vel_sia_observed.metadata(0).set_string('long_name', "x-component of the 'observed' SIA velocities")
450
451 vel_sia_observed.metadata(1).set_name('v_sia_observed')
452 vel_sia_observed.metadata(1).set_string('long_name', "y-component of the 'observed' SIA velocities")
453
454 vel_ssa_observed.copy_from(vel_surface_observed)
455 vel_ssa_observed.add(-1, vel_sia_observed)
456 vecs.add(vel_ssa_observed, writing=True)
457
458 # If the inverse data file has a variable tauc/hardav_true, this is probably
459 # a synthetic inversion. We'll load it now so that it will get written
460 # out, if needed, at the end of the computation in the output file.
461 if PISM.util.fileHasVariable(inv_data_filename, "%s_true" % design_var):
462 design_true = createDesignVec(grid, design_var, '%s_true' % design_var)
463 design_true.regrid(inv_data_filename, True)
464 design_true.read_attributes(inv_data_filename)
465 vecs.add(design_true, writing=saving_inv_data)
466
467 # Establish a logger which will save logging messages to the output file.
468 message_logger = PISM.logging.CaptureLogger(output_filename, 'pismi_log')
469 PISM.logging.add_logger(message_logger)
470 if append_mode or do_restart:
471 message_logger.readOldLog()
472
473 # Prep the output file from the grid so that we can save zeta to it during the runs.
474 if not append_mode:
475 pio = PISM.util.prepare_output(output_filename)
476 pio.close()
477 zeta.write(output_filename)
478
479 # Log the command line to the output file now so that we have a record of
480 # what was attempted
481 PISM.util.writeProvenance(output_filename)
482
483 # Attach various iteration listeners to the solver as needed for:
484
485 # Iteration report.
486 solver.addIterationListener(PISM.invert.ssa.printIteration)
487
488 # Misfit reporting/logging.
489 misfit_logger = PISM.invert.ssa.MisfitLogger()
490 solver.addIterationListener(misfit_logger)
491
492 if inv_method.startswith('tikhonov'):
493 solver.addIterationListener(PISM.invert.ssa.printTikhonovProgress)
494
495 # Saving the current iteration
496 solver.addDesignUpdateListener(PISM.invert.ssa.ZetaSaver(output_filename))
497
498 # Plotting
499 if do_plotting:
500 solver.addIterationListener(InvSSAPlotListener(grid, Vmax))
501 if solver.method == 'ign':
502 solver.addLinearIterationListener(InvSSALinPlotListener(grid, Vmax))
503
504 # Solver is set up. Give the user's prep module a chance to do any final
505 # setup.
506
507 if prep_module is not None:
508 if prep_module.endswith(".py"):
509 prep_module = prep_module[0:-2]
510 exec("import %s as user_prep_module" % prep_module)
511 user_prep_module.prep_solver(solver)
512
513 # Pausing (add this after the user's listeners)
514 if do_pause:
515 solver.addIterationListener(PISM.invert.listener.pauseListener)
516
517 # Run the inverse solver!
518 if do_restart:
519 PISM.logging.logMessage('************** Restarting inversion. ****************\n')
520 else:
521 PISM.logging.logMessage('============== Starting inversion. ==================\n')
522
523 # Try solving
524 reason = solver.solveInverse(zeta_prior, vel_ssa_observed, zeta)
525 if reason.failed():
526 PISM.logging.logError("Inverse solve FAILURE:\n%s\n" % reason.nested_description(1))
527 quit()
528
529 PISM.logging.logMessage("Inverse solve success (%s)!\n" % reason.description())
530
531 (zeta, u) = solver.inverseSolution()
532
533 # It may be that a 'tauc'/'hardav' was read in earlier. We replace it with
534 # our newly generated one.
535 if vecs.has(design_var):
536 design = vecs.get(design_var)
537 design_param.convertToDesignVariable(zeta, design)
538 else:
539 # Convert back from zeta to tauc or hardav
540 design = createDesignVec(grid, design_var)
541 design_param.convertToDesignVariable(zeta, design)
542 vecs.add(design, writing=True)
543
544 vecs.add(zeta, writing=True)
545
546 u.metadata(0).set_name("u_ssa_inv")
547 u.metadata(0).set_string("long_name", "x-component of SSA velocity computed by inversion")
548
549 u.metadata(1).set_name("v_ssa_inv")
550 u.metadata(1).set_string("long_name", "y-component of SSA velocity computed by inversion")
551
552 vecs.add(u, writing=True)
553
554 residual = PISM.model.create2dVelocityVec(grid, name='_inv_ssa_residual')
555 residual.copy_from(u)
556 residual.add(-1, vel_ssa_observed)
557
558 r_mag = PISM.IceModelVec2S()
559 r_mag.create(grid, "inv_ssa_residual", PISM.WITHOUT_GHOSTS, 0)
560
561 r_mag.set_attrs("diagnostic",
562 "magnitude of mismatch between observed surface velocities and their reconstrution by inversion",
563 "m s-1", "m year-1", "inv_ssa_residual", 0)
564 r_mag.metadata().set_number("_FillValue", convert(-0.01, 'm/year', 'm/s'))
565 r_mag.metadata().set_number("valid_min", 0.0)
566
567 r_mag.set_to_magnitude(residual)
568 r_mag.mask_by(vecs.land_ice_thickness)
569
570 vecs.add(residual, writing=True)
571 vecs.add(r_mag, writing=True)
572
573 # Write solution out to netcdf file
574 forward_run.write(output_filename, append=append_mode)
575 # If we're not in append mode, the previous command just nuked
576 # the output file. So we rewrite the siple log.
577 if not append_mode:
578 message_logger.write(output_filename)
579
580 # Save the misfit history
581 misfit_logger.write(output_filename)
582
583
584 if __name__ == "__main__":
585 run()
586
587 # try to stop coverage and save a report:
588 try: # pragma: no cover
589 cov.stop()
590 report = PISM.OptionBool("-report_coverage", "save coverage information and a report")
591 if report:
592 cov.save()
593 cov.html_report(include=["pismi.py"], directory="pismi_coverage")
594 except: # pragma: no cover
595 pass