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